基于OpenCL的射電干涉陣成像網(wǎng)格化算法實(shí)現(xiàn)*

馮 勇， 王 鋒，2，3， 鄧 輝，2， 衛(wèi)守林， 梅 盈，2，3，戴 偉，3， 石聰明

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

中國新一代厘米-分米波射電日像儀——明安圖射電頻譜日像儀(MingantU SpEctral Radioheliograph，MUSER)是基于綜合孔徑成像原理實(shí)現(xiàn)在頻率0.4～15.0 GHz對(duì)太陽進(jìn)行高分辨率(時(shí)間、空間)觀測并成像的射電望遠(yuǎn)鏡[1-2]。

明安圖射電頻譜日像儀主要由天線(天線陣)、接收系統(tǒng)、數(shù)據(jù)處理系統(tǒng)3部分組成。天線陣由低頻陣(MUSER-I，共40面天線)和高頻陣(MUSER-II，共60面天線)兩部分組成，天線陣的拓?fù)浣Y(jié)構(gòu)為螺旋陣[3-4]。在觀測過程中，低頻陣和高頻陣分別每3 ms產(chǎn)生一個(gè)數(shù)據(jù)幀，每個(gè)數(shù)據(jù)幀分別用100 KB和256 KB存儲(chǔ)，數(shù)據(jù)流量分別為32 MB/s和86 MB/s，以每天觀測10小時(shí)為例，每日原始觀測數(shù)據(jù)量分別約為1.2 TB和3.2 TB[5-6]。

如此海量的天文數(shù)據(jù)，對(duì)于數(shù)據(jù)處理工作而言是一個(gè)巨大的挑戰(zhàn)，但挑戰(zhàn)與機(jī)遇并存。在前期研究工作中，為了滿足數(shù)據(jù)處理的需求，研究人員一方面采用分布式計(jì)算技術(shù)，設(shè)計(jì)并實(shí)現(xiàn)了一套分布式數(shù)據(jù)處理執(zhí)行框架OpenCluster[7]；另一方面采用并行計(jì)算技術(shù)并借助高性能計(jì)算設(shè)備，專注于底層算法研究，實(shí)現(xiàn)了一套高效的射電干涉陣數(shù)據(jù)處理管線[8]。特別地，采用消息傳遞接口(Message Passing Interface，MPI)、統(tǒng)一計(jì)算設(shè)備架構(gòu)(Compute Unified Device Architecture，CUDA)以及開放運(yùn)算語言(Open Computing Language，OpenCL)技術(shù)在天文領(lǐng)域有相關(guān)應(yīng)用，并已經(jīng)取得一定的成果[9-11]，證明這些技術(shù)穩(wěn)定性與性能可以滿足天文應(yīng)用軟件開發(fā)的要求。

為了便于明安圖射電頻譜日像儀數(shù)據(jù)處理系統(tǒng)的開發(fā)、部署、應(yīng)用以及推廣，同時(shí)滿足單機(jī)環(huán)境下仍舊能夠進(jìn)行高性能的成像數(shù)據(jù)處理工作，在前期對(duì)OpenCL技術(shù)研究的基礎(chǔ)上，進(jìn)一步采用OpenCL技術(shù)對(duì)射電干涉陣成像中的網(wǎng)格化算法進(jìn)行并行優(yōu)化，在保證網(wǎng)格化算法執(zhí)行效率的同時(shí)，提升算法對(duì)各種硬件平臺(tái)的適應(yīng)性，便于后期數(shù)據(jù)處理軟件的應(yīng)用與推廣。

1 基于OpenCL的網(wǎng)格化算法實(shí)現(xiàn)

1.1 網(wǎng)格化算法

射電干涉陣成像過程是基于綜合孔徑成像原理對(duì)采樣的可見度數(shù)據(jù)進(jìn)行網(wǎng)格化、快速傅里葉變換(Fast Fourier Transformation，F(xiàn)FT)、去卷積(Deconvolution)等操作，最后產(chǎn)生射電源的觀測圖像[12]。

為了借助離散傅里葉變換(Discrete Fourier Transform，DFT)的快速方法——快速傅里葉變換的計(jì)算效率進(jìn)行快速成像，射電干涉陣采集的可見度數(shù)據(jù)必須落在一個(gè)均勻劃分的網(wǎng)格上(如圖1(a))。實(shí)際上，由于UV采樣點(diǎn)分布不均勻，干涉儀得到的是非均勻采樣的可見度數(shù)據(jù)(如圖1(b))，將非均勻采樣的可見度數(shù)據(jù)插值到均勻劃分的網(wǎng)格點(diǎn)上的過程稱為網(wǎng)格化[13-14]。

圖1 均勻劃分的網(wǎng)格(a)與非均勻分布的可見度數(shù)據(jù)示意圖(b)Fig.1 Regular grid(left)and the schematic diagram of non-uniform visibility data(right)

將非均勻采樣可見度數(shù)據(jù)轉(zhuǎn)換成均勻采樣可見度數(shù)據(jù)，即對(duì)未知像素點(diǎn)的數(shù)據(jù)進(jìn)行重采樣，主要有兩類方法：內(nèi)插法和卷積核法。一些內(nèi)插值方法被開發(fā)用于重建非均勻采樣，如最鄰近插值、雙線性插值等。文[15]表明，最優(yōu)的網(wǎng)格化方法是卷積一個(gè)無限擴(kuò)展的sinc函數(shù)，即卷積核法，出于實(shí)際考慮，有限的卷積函數(shù)取代了無限擴(kuò)展的sinc函數(shù)。常用的網(wǎng)格卷積函數(shù)(Gridding Convolution Function，GCF)有截?cái)鄐inc函數(shù)、截?cái)嘀笖?shù)函數(shù)、擬合的球面波函數(shù)。文[16]討論了這些網(wǎng)格卷積函數(shù)在不同尺寸大小下的性能。在日像儀數(shù)據(jù)處理中，選取大小為6的擬合球面波函數(shù)作為網(wǎng)格卷積函數(shù)。卷積網(wǎng)格化并成像過程如圖2。

圖2 卷積網(wǎng)格化并成像Fig.2 Convolutional gridding and imaging

1.2 網(wǎng)格化算法并行化

根據(jù)OpenCL編程規(guī)范，OpenCL程序主要由2部分組成：(1)采用OpenCL語言(類C語言)編寫的內(nèi)核函數(shù)，執(zhí)行在 OpenCL設(shè)備(多核 CPU，GPU，Cell，DSP，F(xiàn)PGA等)上；(2)通過調(diào)用OpenCL的應(yīng)用編程接口(Application Programing Interface，API)管理內(nèi)核函數(shù)的主機(jī)程序，執(zhí)行在主機(jī)(CPU)上。

本文將網(wǎng)格化過程中的關(guān)鍵步驟，由OpenCL語言編寫成內(nèi)核函數(shù)，以便在OpenCL設(shè)備上并行執(zhí)行，并通過在主機(jī)上調(diào)用這些內(nèi)核函數(shù)，從而實(shí)現(xiàn)網(wǎng)格化算法的并行化。主要改寫的內(nèi)核函數(shù)有：gridding＿kernel， correct＿grid＿kernel， normalize＿kernel， point＿symmetric＿kernel， shift＿kernel等。 此外，成像過程的傅里葉變換是通過調(diào)用基于OpenCL實(shí)現(xiàn)的快速傅里葉變換，包含于python的pyfft包中。網(wǎng)格化并行執(zhí)行流程如圖3。

圖3 網(wǎng)格化并行執(zhí)行流程圖Fig.3 The flow diagram of parallel gridding

本文處理的數(shù)據(jù)以二維數(shù)組的形式表示，在采用OpenCL語言編寫內(nèi)核函數(shù)過程中，由OpenCL提供的二維數(shù)組索引get＿global＿id(0)和get＿global＿id(1)，唯一標(biāo)識(shí)二維數(shù)組中的元素，為方便計(jì)算，需借助傳遞給內(nèi)核函數(shù)的二維數(shù)組的參數(shù)(寬度或高度)，將二維數(shù)組按行或按列轉(zhuǎn)換成一維數(shù)組，這樣，二維索引轉(zhuǎn)換成一維索引id，當(dāng)內(nèi)核函數(shù)在并行設(shè)備中執(zhí)行時(shí)，每個(gè)線程處理由索引id唯一標(biāo)識(shí)的二維數(shù)組元素。OpenCL二維索引轉(zhuǎn)一維索引如圖4。

圖4 OpenCL二維索引轉(zhuǎn)一維索引示意圖Fig.4 The schematic diagram of OpenCL two-dimensional index to one-dimensional index

在網(wǎng)格化執(zhí)行過程中，由于可見度數(shù)據(jù)共軛對(duì)稱，在對(duì)可見度數(shù)據(jù)進(jìn)行插值時(shí)，為減少計(jì)算量，先對(duì)一半的可見度數(shù)據(jù)進(jìn)行插值，再通過中心點(diǎn)對(duì)稱方式計(jì)算得到另一半數(shù)據(jù)。將中心點(diǎn)對(duì)稱操作編寫成內(nèi)核函數(shù) point＿symmetric＿kernel， 如圖 5。

圖5 中心點(diǎn)對(duì)稱操作示意圖Fig.5 The schematic diagram of point symmetric operation

point＿symmetric＿kernel偽代碼如下：

在傅里葉變換過程中，為使變換后的圖像是一個(gè)完整的周期，需要進(jìn)行平移操作。將平移操作編寫成內(nèi)核函數(shù)shift＿kernel，如圖6。

圖6 傅里葉變換平移操作示意圖Fig.6 The schematic diagram of Fourier transform shift operation

在網(wǎng)格化執(zhí)行過程中，由于可見度數(shù)據(jù)與一個(gè)網(wǎng)格卷積函數(shù)進(jìn)行了卷積，為消除網(wǎng)格卷積函數(shù)的影響，需進(jìn)行網(wǎng)格校正操作。將網(wǎng)格校正操作編寫成內(nèi)核函數(shù)correct＿grid＿kernel。進(jìn)行網(wǎng)格校正后，需要對(duì)臟圖和臟束進(jìn)行歸一化或標(biāo)準(zhǔn)化，即將臟圖和臟束中所有元素除以臟束中的最大值，從而保證臟束的峰值為1以及臟圖與臟束具有相同量綱，便于后續(xù)進(jìn)行潔化。

1.3 并行加權(quán)

由于UV采樣點(diǎn)太少，對(duì)可見度函數(shù)采樣值進(jìn)行傅里葉逆變換得到的圖像(臟圖)包含一些虛假信息，通過對(duì)UV采樣點(diǎn)賦予不同的權(quán)值改善圖像質(zhì)量的操作稱之為加權(quán)[17]。常用的加權(quán)方式有自然加權(quán)、統(tǒng)一加權(quán)以及Taper，本文將這3種加權(quán)方式改寫成OpenCL內(nèi)核函數(shù)weight＿natural＿kernel，weight＿uniform＿kernel和 weight＿taper＿kernel。

3種加權(quán)方法偽代碼如下：

自然加權(quán)賦予相應(yīng)像素點(diǎn)的權(quán)重為1，統(tǒng)一加權(quán)賦予相應(yīng)像素點(diǎn)的權(quán)重為區(qū)域內(nèi)采樣點(diǎn)頻數(shù)的倒數(shù)，Taper加權(quán)賦予相應(yīng)像素點(diǎn)的權(quán)重為對(duì)應(yīng)高斯函數(shù)值。通過改變像素點(diǎn)的權(quán)重大小，對(duì)可見度數(shù)據(jù)進(jìn)行加權(quán)操作，從而改善圖像質(zhì)量。

2 實(shí)驗(yàn)結(jié)果

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

OpenCL定義不同的公司或廠商為不同的平臺(tái)(Platform)，每個(gè)公司或廠商提供不同的OpenCL設(shè)備(Device)。本文選用的OpenCL設(shè)備有NVIDIA的GPU(GeForce GTX TITAN X，Tesla k20m)和Intel的CPU(Intel Xeon E5-2620 V2)，平臺(tái)及設(shè)備參數(shù)見表1。

表1 平臺(tái)和設(shè)備參數(shù)信息Table 1 The information of platform and device

2.2 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)數(shù)據(jù)來源于明安圖射電頻譜日像儀2015年11月1日12時(shí)8分49秒354毫秒對(duì)太陽的觀測數(shù)據(jù)，原始觀測數(shù)據(jù)經(jīng)過一系列預(yù)處理(壞天線標(biāo)記、星表查詢等)合成為標(biāo)準(zhǔn)天文數(shù)據(jù)格式(UVFITS文件)，數(shù)據(jù)文件被命名為20151101-120849＿354161240.uvfits。通過數(shù)據(jù)處理系統(tǒng)執(zhí)行網(wǎng)格化和快速傅里葉變換等操作，生成相應(yīng)時(shí)刻太陽亮度分布圖像(未經(jīng)過去卷積操作的臟圖，如圖7(a))以及對(duì)應(yīng)的點(diǎn)擴(kuò)散函數(shù)(Point Spread Function，PSF)，點(diǎn)擴(kuò)散函數(shù)也稱為臟束(如圖7(b))。

圖7 MUSER成的臟圖(a)和臟束(b)，觀測時(shí)間是2015年11月1日12：08：49：354，頻率是1.712 5 GHz，極化方式是右旋Fig.7 The dirty image(a)and dirtybeam(b)of MUSER at 12：08：49：354 on November 11th， 2015(Frequency： 1.7125GHz， Polarization： right)

在同一臺(tái)服務(wù)器上，分別在中央處理器(Intel Xeon E5-2620 V2)、圖形處理器(GeForce GTX TITAN X，Tesla k20m)并行設(shè)備上執(zhí)行基于OpenCL實(shí)現(xiàn)的網(wǎng)格化算法和在圖形處理器(GeForce GTX TITAN X，Tesla k20m)上執(zhí)行基于CUDA實(shí)現(xiàn)的網(wǎng)格化算法，并進(jìn)行快速傅里葉變換成圖，分別生成大小為1 024×1 024 pixel和2 048×2 048 pixel的臟圖。網(wǎng)格化并成圖(Gridding＋I(xiàn)FFT)過程執(zhí)行時(shí)間見表2。

從表2可以看出，執(zhí)行基于OpenCL和基于CUDA實(shí)現(xiàn)的網(wǎng)格化算法并成圖，在圖形處理器環(huán)境下消耗的時(shí)間基本相當(dāng)，同時(shí)，基于OpenCL實(shí)現(xiàn)的網(wǎng)格化算法還能在純中央處理器環(huán)境下較快速地執(zhí)行。

表2 Gridding+IFFT過程執(zhí)行時(shí)間表Table 2 The Execution time of Gridding and FFT

2.3 討 論

本文采用并行計(jì)算OpenCL技術(shù)，基于Python語言導(dǎo)入PyOpenCL與PyFFT擴(kuò)展包，實(shí)現(xiàn)了對(duì)日像儀成像過程中的網(wǎng)格化算法以及快速傅里葉變換成圖過程。基于OpenCL實(shí)現(xiàn)的網(wǎng)格化算法具有如下優(yōu)點(diǎn)：

(1)能夠在圖形處理器環(huán)境中執(zhí)行，并且不局限于NVIDIA GPU(實(shí)驗(yàn)條件限制，只選用了NVIDIA GPU和Intel CPU)，執(zhí)行效率與基于CUDA實(shí)現(xiàn)的網(wǎng)格化算法執(zhí)行效率大致相當(dāng)，保證了日像儀成圖的性能；

(2)能夠在多核中央處理器環(huán)境下執(zhí)行，適用于在單機(jī)環(huán)境下進(jìn)行開發(fā)與測試，便于軟件的推廣與應(yīng)用。

綜上所述，采用OpenCL實(shí)現(xiàn)的并行網(wǎng)格化算法，在保證算法執(zhí)行效率的同時(shí)有效地提升了算法對(duì)硬件平臺(tái)的適應(yīng)性。此外，本工作進(jìn)一步驗(yàn)證了OpenCL在天文軟件開發(fā)中的可用性，從CUDA到OpenCL，異構(gòu)系統(tǒng)從NVIDIA GPU＋CPU的異構(gòu)模式轉(zhuǎn)變成并行設(shè)備加中央處理器的異構(gòu)模式。這種異構(gòu)模式的轉(zhuǎn)變將擴(kuò)展并行計(jì)算在天文高性能軟件開發(fā)中的應(yīng)用。

3 下一步工作

本文研究的射電干涉陣成像過程局限于小視場，在大視場成像中w項(xiàng)引起的視場扭曲是不可忽略的，傳統(tǒng)的二維快速傅里葉變換成像將不再適用。因此，下一步將在此研究工作的基礎(chǔ)上，采用OpenCL技術(shù)對(duì)大視場成像中的關(guān)鍵算法，例如w-projection，faceting以及大視場混合算法，進(jìn)行并行優(yōu)化。

致謝：感謝國家天文臺(tái)-阿里云天文大數(shù)據(jù)聯(lián)合研究中心對(duì)本項(xiàng)工作的支持。