梅　盈，鄧　輝，張曉麗，王　鋒1,,4，陳騰達(dá)

(1. 中國科學(xué)院云南天文臺(tái)，云南 昆明　650011；2. 昆明理工大學(xué)，云南 昆明　650500； 3. 廣州大學(xué)天體物理中心，廣東 廣州　510006；4. 中國科學(xué)院大學(xué)，北京　100049)

射電觀測(cè)是研究太陽劇烈活動(dòng)的重要探測(cè)手段。建于中國明安圖的螺旋天線陣列——明安圖射電頻譜日像儀(Mingantu Spectral Radioheliograph, MUSER)在400 MHz～15 GHz的584個(gè)頻點(diǎn)上對(duì)太陽進(jìn)行高分辨率成像，其觀測(cè)將填補(bǔ)日本野邊山日像儀(2個(gè)頻點(diǎn)，17 GHz, 34 GHz)、法國南茜日像儀(5個(gè)頻點(diǎn)，150～450 MHz)、俄羅斯伊爾庫茨克射電日像儀(5.7 GHz)等在厘米、分米波段的觀測(cè)空白，對(duì)太陽活動(dòng)及其對(duì)人類影響的研究起到重要作用[1]。

明安圖射電頻譜日像儀分為低頻陣(MUSER-I, 400 MHz～2 GHz)和高頻陣(MUSER-II, 2 GHz～15 GHz)，其高時(shí)間、高空間和高頻率分辨率給數(shù)據(jù)處理系統(tǒng)提出了較大的挑戰(zhàn)。高頻陣和低頻陣的數(shù)據(jù)接收機(jī)每3 ms分別接收一次包含16通道的觀測(cè)數(shù)據(jù)。在后續(xù)的數(shù)據(jù)處理過程中，如何高效地獲取高質(zhì)量的太陽圖像成為數(shù)據(jù)處理系統(tǒng)的關(guān)鍵問題。

射電干涉陣成圖的主要過程包含加權(quán)(Weighting)、網(wǎng)格化(Gridding)、傅里葉變換、退卷積等主要步驟，其中退卷積是對(duì)臟圖進(jìn)行潔化的過程，是成像過程中最耗時(shí)的部分。對(duì)日像儀成像來說，退卷積過程也是整個(gè)管線設(shè)計(jì)中最令人關(guān)注的問題。本文在細(xì)致調(diào)研相關(guān)工作的基礎(chǔ)上，進(jìn)一步對(duì)射電干涉陣成像中的多尺度退卷積方法進(jìn)行了討論，給出了相應(yīng)的實(shí)現(xiàn)方法。

1　相關(guān)研究工作

射電干涉成像的基本原理是天空實(shí)際亮度是可見度函數(shù)的傅里葉變換。由于在實(shí)際觀測(cè)過程中UV覆蓋的不完全性，采樣函數(shù)是離散的。自1974年H?gbom提出潔化算法的思想后，潔化算法得到了極大的發(fā)展，并在射電望遠(yuǎn)鏡的圖像重建算法中得到廣泛的應(yīng)用。經(jīng)典H?gbom CLEAN算法的關(guān)鍵是為天空實(shí)際亮度提供一個(gè)模型(δ函數(shù))，迭代尋找亮度的峰值位置，并如(1)式進(jìn)行疊加，其中參數(shù)Ip為迭代過程中的峰值亮度，(xp,yp)為峰值亮度的位置，由殘留圖像的峰值亮度確定[2-3]。(2)式中的ID為臟圖亮度，B為臟束，IC為(1)式計(jì)算出的疊加后的亮度值，IR為殘圖的亮度，n為迭代次數(shù)。

(1)

IR(n)=ID-B*IC(n-1).

(2)

H?gbom CLEAN算法的一些變種算法，如Clark, Cotton-Schwab算法、最大熵算法(the Maximum Entropy Method, MEM)隨后被提出。Steer Clean算法被用于日本野邊山日像儀的成像，其潔化效率是H?gbom算法的近10倍[4]。然而，基本的H?gbom CLEAN算法由于在處理點(diǎn)源和展源上表現(xiàn)的優(yōu)勢(shì)，一直被廣泛使用。

但隨著新一代射電干涉儀的不斷發(fā)展，在確?；緷嵒惴ǖ暮唵涡院托旁氡鹊那疤嵯?，多尺度成為提高效率的方法之一。Multi-Resolution CLEAN, Multi-Scale Maximum Entropy, Wavelets CLEAN等多尺度算法有效提高了潔化效率，但存在尺度大小固定和計(jì)算代價(jià)高等問題。

文[3]于2008年提出的多尺度潔化(Multi-Scale CLEAN)算法被認(rèn)為算法直接且具有穩(wěn)定性和收斂性好等特點(diǎn)。多尺度潔化為減少算法的迭代次數(shù)，提高算法的穩(wěn)定性，采用多個(gè)不同尺度大小的卷積核與臟圖進(jìn)行卷積，得到多個(gè)卷積圖像，在多個(gè)圖像中尋找全局亮度峰值。算法的基本流程如下：

(1)選取不同的尺度大小，與臟圖進(jìn)行卷積，得到卷積后的系列圖像；

(2)尋找全局亮度最大點(diǎn)，記錄亮度最大點(diǎn)的位置和所在的尺度大?。?/p>

(3)計(jì)算步驟(2)中的尺度與臟束的卷積，再乘以增益值，存儲(chǔ)計(jì)算結(jié)果及尺度大??；

(4)從(1)中所有的卷積圖像中減去(3)的計(jì)算結(jié)果；

(5)重復(fù)以上4個(gè)步驟直至到達(dá)循環(huán)閾值或達(dá)到預(yù)設(shè)循環(huán)次數(shù)；

(6)將(3)中存儲(chǔ)的結(jié)果與潔束卷積并加上參與圖得到最終的潔圖。

為了控制算法的迭代次數(shù)，提高圖片質(zhì)量，多尺度潔化算法需要指定不同尺度的卷積圖像不同的權(quán)值，每次迭代循環(huán)結(jié)束后，殘圖需要乘以所在尺度的權(quán)值來增加或減少亮度幅值。經(jīng)過系列實(shí)驗(yàn)，(3)式中權(quán)值α=0.6被認(rèn)為是較合適的參數(shù)值[3]。

(3)

其中，S(a)為權(quán)值；a為當(dāng)前尺度大??；a(maxscale)為最大尺度大?。沪翞闄?quán)重參數(shù)。

為了減小迭代算法對(duì)增益的依賴，文[3]選擇拋物線函數(shù)作為多尺度潔化算法的卷積核函數(shù)如下：

m(r,α)=ψ(r)[1-(r/a)2]，

(4)

由于高斯卷積核與(4)式給出的卷積核函數(shù)僅在較高動(dòng)態(tài)范圍下存在較小的差異，一般來說，采用不同尺度的高斯卷積核。以潔束大小為基準(zhǔn)，卷積核大小選取常用的等比數(shù)列，即依次取潔束大小的1、2、4、8倍。

為進(jìn)一步提高處理速度，近幾年利用圖形處理器技術(shù)加快退卷積處理成為主要的方法，如基于圖形處理器的網(wǎng)格化算法使得w-projection的效率較中央處理器提高近百倍*https://arxiv.org/abs/1403.4209。前期基于圖形處理器的H?gbom CLEAN算法較中央處理器下的H?gbom效率提高了十倍以上。但從多尺度潔化算法來看，其并行實(shí)現(xiàn)存在較大的難度。因?yàn)樗阉髯畲笾颠^程雖然可以并行，但需要利用找出的最大值要同時(shí)在所有尺度的圖像中退卷積，這意味著一個(gè)尺度的數(shù)據(jù)處理必須與其它尺度的數(shù)據(jù)處理相交叉，影響了并行的實(shí)現(xiàn)。

2　基于頻域的多尺度潔化

2.1　基本思想

參考文[3]的潔化方法，要加快潔化的效率，關(guān)鍵是能夠同時(shí)進(jìn)行潔化。實(shí)際上，由于H?gbom Clean的基本原理，潔化迭代的過程始終是串行的，要并行的唯一可能性，是同時(shí)能夠獨(dú)立地進(jìn)行多個(gè)潔化迭代。

法國默東天文臺(tái)太陽干涉陣與明安圖射電頻譜日像儀一樣，也是一個(gè)太陽專用的射電干涉陣，其數(shù)據(jù)處理采用了基于頻率帶通濾波的多尺度潔化方法[5]。該方法基于兩個(gè)簡化：(1)沒有過多的展源存在；(2)空間尺度中只有相對(duì)較窄的譜。對(duì)于太陽觀測(cè)來說，完全可以滿足這兩個(gè)簡化需要。該方法與多尺度潔化方法有一定的相似性，但更易于并行實(shí)現(xiàn)。

參考法國默東天文臺(tái)的潔化思想，結(jié)合明安圖射電頻譜日像儀的天線分布和最長基線情況，最終在本文實(shí)現(xiàn)的算法基本原理如下：

(1)在頻域UV平面上，構(gòu)建一系列連續(xù)的帶通濾波器Fk(u,v)，k為濾波器編號(hào)，濾波器的半徑以2的倍數(shù)增長。不同的濾波器對(duì)應(yīng)不同的空間尺度。與法國默東太陽干涉陣不同，在實(shí)現(xiàn)時(shí)每次將計(jì)算最長基線和相應(yīng)的UV分布，反過來計(jì)算相應(yīng)的濾波器帶通頻率。

(2)通過濾波器，將原來的稀疏UV分布D(u,v)進(jìn)一步拆分為若干個(gè)不同尺度的獨(dú)立分布，即：Dk(u,v)=Fk(u,v)×D(u,v)；

(3)對(duì)于一個(gè)不同尺度的UV分布，通過逆傅里葉變換得到臟圖：Ik=FT[Dk(u,v)]，然后進(jìn)行常規(guī)的潔化處理；

(4)最后的潔圖等于各尺度潔圖的累加。

2.2　尺度的考慮

根據(jù)5個(gè)不同尺度的大小，將臟圖分成不同的頻率段，隨著k值增大，尺度依次減小，直到覆蓋所有的UV值。當(dāng)k=0時(shí)，濾波器為高斯濾波。其尺度為最長UV的20%，從成像來看，這也是最大的尺度。k=1、2、3、4時(shí)，濾波器在頻率域呈現(xiàn)一個(gè)環(huán)狀，實(shí)現(xiàn)一個(gè)帶通濾波。無論分成多少級(jí)，但本質(zhì)上要確保這些濾波器和所合成的濾波器在頻率上連續(xù)，覆蓋全部UV點(diǎn)，如圖1和圖2中的最后一張圖。在使用高斯濾波的情況下，最終合成的濾波器在邊緣逐漸減小的濾波性質(zhì)，正好減小了由于較大的UV處的稀疏導(dǎo)致的成圖不可靠信息。

3　算法的圖形處理器實(shí)現(xiàn)及效率

3.1　圖形處理器實(shí)現(xiàn)

明安圖射電頻譜日像儀成圖的圖形處理器實(shí)現(xiàn)采用GPU-CUDA(Compute Unified Device Architecture)架構(gòu)，使用Python編程語言，PyCUDA，Scikit-CUDA作為CUDA并行編程接口庫。實(shí)現(xiàn)了數(shù)據(jù)處理的單機(jī)命令行模式及分布式運(yùn)行模式[6]，其中數(shù)據(jù)預(yù)處理(數(shù)據(jù)格式轉(zhuǎn)換[7]、相位校準(zhǔn)、異常數(shù)據(jù)標(biāo)記等)在中央處理器環(huán)境下運(yùn)行，成圖及潔化在圖形處理器環(huán)境下運(yùn)行，包含8臺(tái)服務(wù)器的高性能計(jì)算機(jī)群(CPUs + GPUs)已在觀測(cè)站搭建好。成像流程如圖3，其中日像儀多尺度潔化算法描述如表1。以低頻陣觀測(cè)數(shù)據(jù) “2015-11-01 04:08:49. 354161240(UTC), 1.712 5 GHz，右旋” 為例，最終潔化圖像如圖4(a)。在實(shí)驗(yàn)過程中，常用天文學(xué)軟件應(yīng)用包(the Common Astronomy Software Applications package, CASA)用于驗(yàn)證數(shù)據(jù)處理各個(gè)階段的正確性。將數(shù)據(jù)處理系統(tǒng)生成的UVFITS文件通過CASA轉(zhuǎn)換成MS文件并使用多尺度潔化成圖可得到基本一致的成圖結(jié)果。為說明當(dāng)前處理的正確性，文中給出了日本野邊山天文臺(tái)已經(jīng)公布的觀測(cè)處理結(jié)果圖4(b)作為對(duì)比。由于當(dāng)前的數(shù)據(jù)處理過程中沒有進(jìn)行日軸方位角改正，圖像中存在一些誤差，在后期的工作中將進(jìn)行改正。

圖1濾波器設(shè)計(jì)
Fig.1Filters profiles alongu(-) andv(-) axes

圖2UV濾波器設(shè)計(jì)
Fig.2Filters in the UV-plane

圖3基于圖形處理器的成像流程
Fig.3GPU-based implementation of imaging

表1　明安圖射電頻譜日像儀多尺度潔化Table 1　Multi-Scale CLEAN for MSUER

3.2　算法效率

由于由多個(gè)濾波器得到的圖像同時(shí)進(jìn)行潔化，每個(gè)濾波后的圖像有不同的迭代次數(shù)，且潔化算法的迭代次數(shù)動(dòng)態(tài)改變，因此本文僅以具體實(shí)例說明多尺度帶通濾波的性能。在NVIDIA Corporation GM200測(cè)試環(huán)境下，對(duì)于處理圖4中的2015-11-01 1.712 5 GHz的一張1 024 × 1 024的圖像，使用基于帶通濾波的多尺度潔化時(shí)間為1.424 s，前期實(shí)驗(yàn)中的H?gbom CLEAN算法[8]的潔化時(shí)間為4.037 48 s。經(jīng)過多次潔化處理實(shí)驗(yàn)，表明當(dāng)前的多尺度潔化能有效將成圖效率提高近3倍，同時(shí)可以減少迭代次數(shù)，對(duì)數(shù)據(jù)處理流水線性能的提升有重要意義。

圖4(a) 明安圖射電頻譜日像儀潔化圖(Multi-scale)(MUSER-I at 2015-11-01 04:08:49. 354161240 (UTC)，頻率: 1.712 5 GHz，右旋)； (b) 野邊山天文臺(tái)潔化圖(2015-11-01 at 04: 10: 00 (UTC), 頻率: 17 GHz (R + L)

Fig.4(a) Clean image of MUSER on 2015-11-01 at 04:08:49. 354161240 (UTC), frequency: 1.7125GHz, polarization: right); (b) Clean image from Nobeyama on 2015-11-01 at 04:10:00 (UTC), frequency: 17GHz (R + L)

4　總　結(jié)

本文在分析多尺度潔化算法原理的基礎(chǔ)上，實(shí)現(xiàn)了基于帶通濾波的多尺度潔化算法，較文[3]的多尺度潔化算法在全局不同尺度的圖像中的迭代而言，完全實(shí)現(xiàn)了不同尺度的并行。實(shí)驗(yàn)選取了針對(duì)明安圖射電頻譜日像儀成像的濾波器、尺度大小及不同尺度的取值。同時(shí)，基于圖形處理器的實(shí)現(xiàn)表明，基于帶通濾波的多尺度潔化算法比同等條件下的H?gbom CLEAN性能提高了近3倍，并且可有效減少迭代次數(shù)。本文的研究結(jié)果將有效提高整個(gè)數(shù)據(jù)處理系統(tǒng)的性能，并為射電干涉成像提供了有價(jià)值的參考。