寬視場(chǎng)成像網(wǎng)格化算法中w-plane最優(yōu)經(jīng)驗(yàn)值研究*

于曉雨，鄧 輝,2，梅 盈，衛(wèi)守林，石聰明，王 威，戴 偉，王 鋒,2

(1. 昆明理工大學(xué)云南省計(jì)算機(jī)技術(shù)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650500; 2. 廣州大學(xué)天體物理中心/物理與電子工程學(xué)院，廣東 廣州 510006；3. 中國(guó)科學(xué)院國(guó)家天文臺(tái)，北京 100101)

平方千米陣列(Square Kilometer Array, SKA)是目前國(guó)際上建造的最大綜合孔徑射電望遠(yuǎn)鏡[1]，它由數(shù)千個(gè)反射面天線和多達(dá)一百萬(wàn)個(gè)低頻天線組成，總接收面積達(dá)到一平方千米。

平方千米陣列科學(xué)數(shù)據(jù)處理(Science Data Processor, SDP)是望遠(yuǎn)鏡設(shè)計(jì)中的一個(gè)重要組成部分?？茖W(xué)數(shù)據(jù)處理專(zhuān)注于設(shè)計(jì)硬件平臺(tái)、軟件和算法處理射電望遠(yuǎn)鏡的觀測(cè)數(shù)據(jù)[2]。科學(xué)數(shù)據(jù)處理的設(shè)計(jì)面臨的挑戰(zhàn)是多方面的，需要處理來(lái)自望遠(yuǎn)鏡的海量觀測(cè)數(shù)據(jù)。SKA-1(約占10%的總體規(guī)模)數(shù)據(jù)量將達(dá)到300 PB/年。為了滿(mǎn)足平方千米陣列科學(xué)研究開(kāi)展的需要，科學(xué)數(shù)據(jù)處理在數(shù)據(jù)處理管線設(shè)計(jì)中，需要充分考慮數(shù)據(jù)處理的性能與質(zhì)量。由于科學(xué)數(shù)據(jù)處理需要的計(jì)算力高出現(xiàn)有最大的天文數(shù)據(jù)處理系統(tǒng)兩個(gè)數(shù)量級(jí)，數(shù)據(jù)處理關(guān)鍵算法的運(yùn)行速率要達(dá)到百萬(wàn)兆級(jí)的運(yùn)算水平。SKA-2的計(jì)算力需求估計(jì)達(dá)到4 exaFLOPS的水平。因此，開(kāi)展對(duì)科學(xué)數(shù)據(jù)處理中關(guān)鍵算法的研究，提高數(shù)據(jù)處理的性能具有顯著的意義。

網(wǎng)格化(Gridding)算法是綜合孔徑望遠(yuǎn)鏡數(shù)據(jù)處理中的一種重要算法，對(duì)計(jì)算的消耗占據(jù)了很大比重。對(duì)網(wǎng)格化算法進(jìn)行性能優(yōu)化，對(duì)科學(xué)數(shù)據(jù)處理有顯著的價(jià)值和意義。低頻射電望遠(yuǎn)鏡大視場(chǎng)成像受到很多條件的制約，其中最重要的一個(gè)影響因素是非共面基線效應(yīng)。在大視場(chǎng)與非共面基線效應(yīng)下，對(duì)可見(jiàn)度數(shù)據(jù)直接進(jìn)行傅里葉變換會(huì)導(dǎo)致圖像嚴(yán)重畸變，因此，在大視場(chǎng)成圖時(shí)，解決w項(xiàng)的問(wèn)題至關(guān)重要。目前，應(yīng)對(duì)大視場(chǎng)與非共面基線效應(yīng)的算法有：faceting[3]，三維傅里葉變換[4]，w-projection[5]，w-stacking[6]和warped snapshots[4]等算法以及其他復(fù)合算法。在這些方法中，w-projection和w-stacking是目前應(yīng)用相對(duì)較廣的網(wǎng)格化算法。其中，w-projection算法在犧牲較大內(nèi)存的情況下，具有優(yōu)越的計(jì)算速度和誤差控制。該方法不僅可以使用傳統(tǒng)的潔化方法去卷積，還可以選擇多尺度潔化方法，在高動(dòng)態(tài)范圍下，有較大的優(yōu)勢(shì)，相位中心基本不失真。w-stacking算法在圖像空間以消耗更多的內(nèi)存為代價(jià)處理傅里葉空間中w項(xiàng)，計(jì)算速度優(yōu)于w-projection[7]。

在w-projection和w-stacking算法的具體實(shí)現(xiàn)過(guò)程中，w-plane是一個(gè)重要參數(shù)。在w的最大值和最小值之間分出多個(gè)w-plane，對(duì)不同w-plane的可見(jiàn)度數(shù)據(jù)使用不同的卷積核進(jìn)行卷積運(yùn)算。因此，w-plane的數(shù)量直接影響算法的計(jì)算速度、內(nèi)存開(kāi)銷(xiāo)以及成圖質(zhì)量。本文主要對(duì)w-projection和w-stacking算法在實(shí)現(xiàn)與實(shí)際運(yùn)行中的w-plane進(jìn)行實(shí)驗(yàn)研究，利用SKA-1低頻陣臺(tái)站數(shù)據(jù)和澳大利亞SKA探路者(Australian SKA Pathfinder, ASKAP)軟件包進(jìn)行模擬觀測(cè)，以期掌握w-plane對(duì)整體計(jì)算性能的影響，為后期的平方千米陣列科學(xué)數(shù)據(jù)處理建設(shè)以及國(guó)內(nèi)科學(xué)數(shù)據(jù)中心建設(shè)中的管線系統(tǒng)設(shè)計(jì)提供參考。

1 w-plane

低頻射電望遠(yuǎn)鏡大視場(chǎng)成像最重要的一個(gè)影響因素是非共面基線效應(yīng)，在非共面基線效應(yīng)的影響下，w的值遠(yuǎn)大于1，因此直接對(duì)可見(jiàn)度數(shù)據(jù)進(jìn)行二維傅里葉變換的條件已經(jīng)不滿(mǎn)足，針對(duì)這一問(wèn)題，天文學(xué)家提出了一系列算法用來(lái)應(yīng)對(duì)非共面基線效應(yīng)造成的圖像畸變。w-projection算法通過(guò)將不同w值的可見(jiàn)度數(shù)據(jù)投影到w=0的平面，從而消除w項(xiàng)的影響，有效提高成圖質(zhì)量。w-stacking算法能夠在像平面矯正非共面基線效應(yīng)引起的相位偏移，從而有效還原天空亮度分布。

在射電綜合孔徑成像處理過(guò)程中，可見(jiàn)度數(shù)據(jù)與天空亮度之間的關(guān)系：

(1)

w-projection算法將不同w值的可見(jiàn)度數(shù)據(jù)投影到w=0的平面，(1)式簡(jiǎn)化為

(2)

經(jīng)過(guò)簡(jiǎn)單的指數(shù)運(yùn)算，(1)式變?yōu)?3)式和(4)式：

(3)

(4)

(3)式、(4)式可以改寫(xiě)為(5)式，其中卷積核通過(guò)(6)式近似計(jì)算得到：

V(u,v,w)=G^(u,v,w)*V(u,v,w=0),

(5)

(6)

以上各式中，V(u,v,w)是w非零的可見(jiàn)度數(shù)據(jù)；I(l,m)表示天空亮度分布；V(u,v,w=0)表示

w=0的可見(jiàn)度數(shù)據(jù)；G(l,m,n)表示w項(xiàng)；G^(l,m,n)表示卷積核。

上述過(guò)程顯示，在w-projection算法中，卷積核是算法的核心。理想情況下，位于不同w-plane的可見(jiàn)度數(shù)據(jù)，在計(jì)算過(guò)程中應(yīng)該采用不同的卷積核。在實(shí)際計(jì)算過(guò)程中，一般計(jì)算w的最大和最小值，在最大和最小值之間平均分成多個(gè)w-plane，然后對(duì)每個(gè)w-plane按照(6)式計(jì)算卷積核并保存，最后對(duì)每個(gè)可見(jiàn)度數(shù)據(jù)與其距離最近的w-plane的卷積核進(jìn)行卷積，從而實(shí)現(xiàn)到w=0平面的投影。

顯而易見(jiàn)，w-plane的值越大，投影效果越好，但是計(jì)算量和內(nèi)存消耗大；w-plane的值越小，計(jì)算量和內(nèi)存消耗小，但是投影效果欠佳。當(dāng)w-plane的數(shù)量為0，則直接退化為二維傅里葉變換。文[5]在前期研究中取最大的w值的平方根作為經(jīng)驗(yàn)值。

2 w-stacking 算法中w-plane

w-stacking算法的提出是為了消除非共面基線效應(yīng)對(duì)大視場(chǎng)成像的影響，但是不同于w-projection算法，w-stacking算法并沒(méi)有消除w項(xiàng)，而是對(duì)其進(jìn)行了修正。

根據(jù)(1)式將uv空間轉(zhuǎn)換到lm空間得到(7)式：

(7)

對(duì)(7)式兩邊在w范圍內(nèi)進(jìn)行積分得到(8)式：

(8)

其中，wmax和wmin分別為w的最大值和最小值。(8)式等號(hào)右側(cè)包括兩部分：(1)二維傅里葉變換；(2)對(duì)w的復(fù)數(shù)積分。對(duì)(8)式等號(hào)右側(cè)離散化得到(9)式：

(9)

w-stacking算法的實(shí)現(xiàn)步驟如下：

(1)根據(jù)可見(jiàn)度數(shù)據(jù)w的范圍設(shè)定w-plane的細(xì)分粒度；

(2)依據(jù)每個(gè)可見(jiàn)度數(shù)據(jù)的w值，將可見(jiàn)度數(shù)據(jù)網(wǎng)格化到某一w-plane；

(3)對(duì)每個(gè)w-plane的可見(jiàn)度數(shù)據(jù)進(jìn)行二維傅里葉變換得到相應(yīng)圖像；

(4)根據(jù)w值對(duì)各個(gè)w-plane的圖像進(jìn)行相位矯正；

(5)加權(quán)疊加所有圖像，得到最終圖像。

在w-stacking算法實(shí)現(xiàn)過(guò)程中，w-plane分的過(guò)粗導(dǎo)致產(chǎn)生偽影，分的過(guò)細(xì)又增大計(jì)算量和內(nèi)存消耗。同時(shí)，w-stacking算法將不同w-plane上的可見(jiàn)度數(shù)據(jù)看成是相互獨(dú)立的，各個(gè)面的數(shù)據(jù)可以并行進(jìn)行二維傅里葉變換和相位改正，因此，w-stacking算法在并行實(shí)現(xiàn)方面具有優(yōu)勢(shì)。

3 觀測(cè)模擬仿真實(shí)驗(yàn)

3.1 實(shí)驗(yàn)環(huán)境

為了驗(yàn)證w-plane的影響，獲得較為理想的w-plane值，在確保成像質(zhì)量的情況下降低成像處理時(shí)間，本文進(jìn)行了一系列實(shí)驗(yàn)。數(shù)據(jù)處理軟件采用澳大利亞SKA探路者的數(shù)據(jù)處理軟件系統(tǒng)[8]。澳大利亞SKA探路者通過(guò)名為ASKAPsoft的軟件系統(tǒng)進(jìn)行望遠(yuǎn)鏡觀測(cè)控制和數(shù)據(jù)處理[9]，ASKAPsoft框架由3部分組成：望遠(yuǎn)鏡觀測(cè)系統(tǒng)、處理中心和科學(xué)數(shù)據(jù)存檔。澳大利亞SKA探路者通過(guò)bash shell腳本組合ASKAPsoft中的程序形成數(shù)據(jù)處理管線，主要的數(shù)據(jù)處理管線有數(shù)據(jù)攝取管線、定標(biāo)管線和成圖管線[10]。

本實(shí)驗(yàn)在8臺(tái)漢柏服務(wù)器上完成，具體硬件配置如表1。

3.2 模擬觀測(cè)與成像處理

為了分析w-plane對(duì)w-projection算法和w-stacking算法的影響，采用模擬觀測(cè)、成像、對(duì)后續(xù)圖像進(jìn)行分析的方法驗(yàn)證不同的w-plane造成的影響。

表1 硬件設(shè)備參數(shù)Table 1 The parameters of device

3.2.1 模擬觀測(cè)

實(shí)驗(yàn)采用SKA1-Low望遠(yuǎn)鏡的512個(gè)天線臺(tái)站坐標(biāo)數(shù)據(jù)，天線布局如圖1。天空模型由73個(gè)點(diǎn)源構(gòu)成，整體形狀為 “米” 字，如圖2(a)。實(shí)驗(yàn)采用16個(gè)計(jì)算節(jié)點(diǎn)，每個(gè)計(jì)算節(jié)點(diǎn)處理一個(gè)通道的數(shù)據(jù)，通道帶寬1 MHz，因此觀測(cè)頻率覆蓋范圍16 MHz，積分時(shí)間間隔150 s，模擬觀測(cè)時(shí)長(zhǎng)6 h。SKA1-Low天線頻率覆蓋范圍50～350 MHz，本實(shí)驗(yàn)中心參考頻率為100 MHz，設(shè)置16個(gè)帶寬為1 MHz的通道，因此觀測(cè)頻率范圍為92～108 MHz。

使用模擬觀測(cè)程序?qū)μ炜漳Ｐ瓦M(jìn)行16個(gè)通道的模擬觀測(cè)，得到18 837 504行5.5 GB的模擬可見(jiàn)度數(shù)據(jù)。

圖1 天線分布
Fig.1 Antennas layout

3.2.2 成像處理

使用w-projection算法和w-stacking算法，對(duì)模擬觀測(cè)獲得的可見(jiàn)度數(shù)據(jù)16個(gè)通道進(jìn)行連續(xù)譜成圖。ASKAP成圖管線在8個(gè)計(jì)算節(jié)點(diǎn)上并行運(yùn)算，每個(gè)計(jì)算節(jié)點(diǎn)處理兩個(gè)通道。實(shí)驗(yàn)過(guò)程中不斷改變w-plane的數(shù)量，得到不同w-plane設(shè)置下的圖像，并記錄每次成圖時(shí)間供后續(xù)分析使用。

圖2(b)是二維傅里葉處理得到的圖像，可以看到w項(xiàng)影響導(dǎo)致的拖尾效應(yīng)(Smearing)非常明顯。圖3是w-plane數(shù)量設(shè)置為171的情況下，w-projection算法和w-stacking算法處理得到的圖像，拖尾效應(yīng)得到抑制，但邊緣的模擬星并沒(méi)有被很好地還原。

3.2.3 圖像質(zhì)量評(píng)價(jià)方法

為了評(píng)價(jià)不同w-plane數(shù)量下的兩種算法的成像質(zhì)量，研究采用了SELAVY軟件包，該軟件包是Duchamp星源搜索算法在ASKAP軟件中的實(shí)現(xiàn)模塊，用來(lái)搜索觀測(cè)結(jié)果中的星數(shù)量。由于在仿真觀測(cè)時(shí)已經(jīng)預(yù)設(shè)了總的星數(shù)為73，同時(shí)所有星的位置已知。因此，通過(guò)SELAVY搜索出的星與已知天空模型星進(jìn)行對(duì)比，利用正確還原的星數(shù)量來(lái)快速判斷成像的質(zhì)量。圖4(a)顯示了隨著w-plane數(shù)量的增長(zhǎng)，兩種算法處理時(shí)間的變化。圖4(b)顯示了隨著w-plane數(shù)量的增長(zhǎng)，SELAVY搜索出點(diǎn)源數(shù)量的變化。

圖2(a) 天空模型圖; (b) 二維傅里葉變換成像結(jié)果
Fig.2(a) Sky model image; (b) The result of two-dimensional Fourier transform imaging

圖3(a) w-projection算法的成圖結(jié)果； (b) w-stacking算法的成圖結(jié)果
Fig.3(a) The result of the w-projection imaging; (b) The result of the w-stack imaging

圖4(a) 計(jì)算時(shí)間與w-plane的關(guān)系； (b) SELAVY搜索出點(diǎn)源數(shù)量與w-plane的關(guān)系

Fig.4(a) The relation between the running time and w-plane;(b) The relation between the number of searched sources and w-plane

4 結(jié)果分析

本實(shí)驗(yàn)最大基線長(zhǎng)度81 602 m，最大w約為28 000個(gè)波長(zhǎng)。根據(jù)文[5]研究的結(jié)果，w-plane值一般設(shè)置為最大w的平方根，也就是167。

SELAVY軟件包對(duì)兩種算法在不同的w-plane時(shí)成圖結(jié)果進(jìn)行星像統(tǒng)計(jì)。從圖4(b)可以看到在w-plane較小的情況下，w-stacking算法成圖得到星的數(shù)量多于w-projection算法的結(jié)果，但是隨著w-plane的增加兩者的成圖質(zhì)量近似相同。

5 結(jié) 論

綜上所述，針對(duì)w-projection算法和w-stacking算法在管線中的應(yīng)用，有如下結(jié)論：

本文的工作表明，平方千米陣列數(shù)據(jù)處理管線應(yīng)根據(jù)不同的需求選擇不同的網(wǎng)格化算法，后續(xù)工作以提高運(yùn)算速度、提高算法的并行性為研究重點(diǎn)，以獲得更好的處理能力。