MUSER可見(jiàn)度數(shù)據(jù)積分方法與實(shí)現(xiàn)*

賴(lài) 鋮，梅 盈,2，鄧 輝，王 鋒,2，戴 偉

(1. 昆明理工大學(xué)云南省計(jì)算機(jī)技術(shù)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650500；2. 中國(guó)科學(xué)院云南天文臺(tái)，云南 昆明 650011)

射電天文是天文學(xué)的一個(gè)分支，它使用射電望遠(yuǎn)鏡系統(tǒng)在無(wú)線電波段研究來(lái)自深空(包括各類(lèi)天體)的射電波[1]。為推動(dòng)射電觀測(cè)技術(shù)的發(fā)展，填補(bǔ)國(guó)際上對(duì)太陽(yáng)耀斑能量初始釋放區(qū)分米波段高分辨射電成像觀測(cè)的空白[2-3]，我國(guó)已建成在0.4～15 GHz范圍內(nèi)的厘米-分米波日像儀——明安圖射電頻譜日像儀。它可以在超寬頻帶下同時(shí)以高時(shí)間、高空間和高頻率分辨率觀測(cè)太陽(yáng)的動(dòng)力學(xué)性質(zhì)，探索太陽(yáng)劇烈活動(dòng)的起源[4]。

MUSER天線陣包含低頻陣(MUSER-I)和高頻陣(MUSER-II)。低頻陣共有40面4.5 m天線，觀測(cè)頻率為400 MHz～2 GHz。低頻陣的觀測(cè)分為4個(gè)波段，兩種極化方式(左旋、右旋)，其中每個(gè)波段包含16個(gè)通道。高頻陣共60面2 m天線，觀測(cè)頻率范圍為2 GHz～15 GHz，在33個(gè)波段和兩種極化方式上觀測(cè)，每個(gè)波段包含16個(gè)通道[5]。高頻陣和低頻陣每3 ms產(chǎn)生一幀觀測(cè)數(shù)據(jù)(一個(gè)波段、一個(gè)極化方向上16個(gè)通道的數(shù)據(jù))。如何將原始觀測(cè)的可見(jiàn)度數(shù)據(jù)經(jīng)過(guò)一系列校準(zhǔn)和計(jì)算得到用于科學(xué)研究的太陽(yáng)圖像，是數(shù)據(jù)處理的關(guān)鍵。

目前，國(guó)外用于太陽(yáng)射電成像觀測(cè)的設(shè)備主要有日本野邊山的太陽(yáng)射電日像儀，積分時(shí)間為25 ms[6]；法國(guó)南茜米波太陽(yáng)射電日像儀，積分時(shí)間1 s①；西伯利亞太陽(yáng)射電日像儀，積分時(shí)間0.336 s[7]。這些射電日像儀積分處理過(guò)程一般直接在接收機(jī)完成，因此時(shí)間分辨率較低。我國(guó)位于內(nèi)蒙古明安圖觀測(cè)站的明安圖射電頻譜日像儀已經(jīng)建成，滿足實(shí)時(shí)和事后科學(xué)數(shù)據(jù)處理的軟件已開(kāi)發(fā)完成[8-9]。為進(jìn)一步提高成圖質(zhì)量，本文通過(guò)對(duì)成圖基本原理的分析，進(jìn)一步提出基于UV平均積分和UV覆蓋積分兩種方法提高太陽(yáng)圖像的信噪比。

1 成像原理

1.1 綜合孔徑成像原理

相關(guān)干涉儀是將間隔一定距離的兩面天線接收的天體輻射分別傳送到相關(guān)器中進(jìn)行相關(guān)處理。由干涉成像原理可知，干涉儀的輸出響應(yīng)和天空亮度存在傅里葉變換關(guān)系，假設(shè)天空亮度分布為I(l,m)，干涉儀輸出為V(u,v)，則有如下關(guān)系[10]：

V(u,v)=∑∑I(l,m)ei2π(ul+vm).

(1)

將離散求和轉(zhuǎn)換成連續(xù)積分可以得出如下公式[10]：

(2)

其中，u,v為投影平面也稱(chēng)作空間頻率；l,m為平時(shí)的天空面。 (2)式表明干涉儀的輸出是亮度分布傅里葉變換的一個(gè)(u,v)分量，隨著u,v的不同可以得到不同空間頻率分量的值，所有的空間頻率分量總稱(chēng)復(fù)可見(jiàn)度函數(shù)。對(duì)復(fù)可見(jiàn)度函數(shù)進(jìn)行傅里葉逆變換得到天空亮度分布，即圖像，這被稱(chēng)為綜合孔徑成像原理[10]。

(3)

1.2 MUSER積分的基本原理

當(dāng)前數(shù)據(jù)處理系統(tǒng)流水線已完成從原始觀測(cè)數(shù)據(jù)格式到FITS文件格式的轉(zhuǎn)換，觀測(cè)數(shù)據(jù)校準(zhǔn)到臟圖以及潔化的成圖過(guò)程[8]。以高頻陣在2016年7月5日4時(shí)4分38秒的觀測(cè)數(shù)據(jù)為例，生成臟圖(圖1(a))和潔化圖(圖1(b))。

圖1 MUSER原始臟圖(a)和潔化圖(b)
Fig.1 Dirty and clean image without integration (MUSER-II, 2016-07-05 04:04:38 (UTC))

可以看出，對(duì)單幀采樣數(shù)據(jù)成像后的太陽(yáng)圖像信噪比很低，不能清楚地看出太陽(yáng)爆發(fā)的位置及其輪廓。臟圖潔化后得到的潔圖，同樣也不能清楚地分析出爆發(fā)位置。為提高日像儀重建的太陽(yáng)圖像的信噪比，得到質(zhì)量較好的成像結(jié)果，需要對(duì)采樣數(shù)據(jù)進(jìn)行積分處理，改善圖像質(zhì)量，提高信噪比[12]。本文根據(jù)科學(xué)研究的需要，以及積分時(shí)段要求，從兩方面分析，研究了UV平均積分(可見(jiàn)度數(shù)據(jù)求和再求平均)、UV覆蓋積分(利用地球自轉(zhuǎn)增加UV覆蓋)兩種方法。

1.2.1 基于UV平均的積分

明安圖射電頻譜日像儀選擇觀測(cè)模式為快照模式，積分時(shí)間為3 ms進(jìn)行觀測(cè)。如此高的時(shí)間分辨率情況下，為了提高空間分辨率，達(dá)到改善圖像質(zhì)量的效果，采取對(duì)3 ms觀測(cè)得到的UV及其輸出數(shù)據(jù)平均的方法實(shí)現(xiàn)。

以高頻陣為例，其中單幀積分采樣時(shí)間tc，假設(shè)對(duì)UV積分時(shí)間為tI，單幀采樣函數(shù)S(u,v)，采樣輸出為V(u,v)，則UV積分表示如下：

(4)

(4)式做法即將采樣的復(fù)可見(jiàn)函數(shù)進(jìn)行平均，同時(shí)采樣函數(shù)也進(jìn)行平均，用平均采樣函數(shù)和復(fù)可見(jiàn)度函數(shù)的方法實(shí)現(xiàn)UV平均積分。

1.2.2 基于UV覆蓋增加的積分

UV覆蓋對(duì)成圖起至關(guān)重要的作用，由于采樣函數(shù)S(u,v)的傅里葉變換即得到臟束，也稱(chēng)綜合束，而臟束又直接影響成圖，因此可通過(guò)增加UV覆蓋提高圖像質(zhì)量。增加UV 覆蓋有3種方式，(1)增加陣列望遠(yuǎn)鏡個(gè)數(shù)；(2)增加子孔徑有效半徑；(3)利用地球自轉(zhuǎn)[12]。增加陣列望遠(yuǎn)鏡個(gè)數(shù)對(duì)成本有直接需求，在成本固定的情況下不考慮，增加子孔徑有效半徑可能導(dǎo)致采樣率不足。在現(xiàn)有成本下通過(guò)地球自轉(zhuǎn)增加UV覆蓋是一種非常直接可行的方法。

隨著地球自轉(zhuǎn)UV發(fā)生改變，而短時(shí)間內(nèi)UV的變化很小，因此假設(shè)每間隔tinterval增加UV覆蓋，總積分時(shí)間ta，則增加UV覆蓋后的采樣函數(shù)和采樣輸出表示為

(5)

上述兩種方式均可以提高成圖的質(zhì)量，將兩種方法混合到一起同樣可以提升圖像質(zhì)量，由此提出先進(jìn)行UV平均積分再增加UV覆蓋積分的方法。

1.2.3 基于UV平均和增加覆蓋的混合積分

假設(shè)UV積分時(shí)間為tI，可以知道在該段時(shí)間內(nèi)的采樣函數(shù)為StI(ui,vi)，采樣輸出VtI(ui,vi)，每間隔tinterval增加UV覆蓋，總積分時(shí)間ta，混合積分采樣函數(shù)及其采樣輸出可以表示為

(6)

2 積分過(guò)程實(shí)現(xiàn)

明安圖射電頻譜日像儀的成圖過(guò)程：目標(biāo)源放射出天體輻射，然后干涉儀對(duì)其進(jìn)行采樣，采樣得到的原始數(shù)據(jù)經(jīng)過(guò)預(yù)處理相關(guān)校準(zhǔn)[13]、對(duì)應(yīng)極化和通道，之后進(jìn)行積分處理，再網(wǎng)格化傅里葉逆變換得到臟圖，最后經(jīng)過(guò)潔化算法處理得到潔圖，過(guò)程如圖2。

圖2 MUSER積分成圖過(guò)程
Fig.2 The procedure of MUSER integration

2.1 積分操作流程

明安圖射電頻譜日像儀的每3 ms接收一幀觀測(cè)數(shù)據(jù)，對(duì)這些數(shù)據(jù)進(jìn)行積分操作得到積分?jǐn)?shù)據(jù)，積分?jǐn)?shù)據(jù)是觀測(cè)得到并且經(jīng)過(guò)預(yù)處理相關(guān)校準(zhǔn)、對(duì)應(yīng)極化和通道的可見(jiàn)度數(shù)據(jù)。UV平均積分具體操作流程如圖3。

圖3 可見(jiàn)度數(shù)據(jù)平均的計(jì)算
Fig.3 Calculation of the average of visibilities

高頻陣中每一幀完整的數(shù)據(jù)包含33小幀。UV積分對(duì)連續(xù)的每一幀數(shù)據(jù)中UU，VV以及采樣得到的實(shí)部、虛部數(shù)據(jù)進(jìn)行累加再求平均，得到的UV平均積分?jǐn)?shù)據(jù)量和單幀大小數(shù)量一樣多，例如在高頻陣中采樣函數(shù)UV點(diǎn)為3 540個(gè)采樣得到的實(shí)部虛部3 540個(gè)，則經(jīng)過(guò)UV平均積分得到的數(shù)據(jù)同樣是3 540個(gè)UV點(diǎn)及其采樣值。

實(shí)現(xiàn)UV覆蓋積分具體操作流程如圖4。對(duì)高頻陣中數(shù)據(jù)進(jìn)行UV覆蓋積分同樣用完整的數(shù)據(jù)幀積分，不同的是由于短時(shí)間內(nèi)采樣函數(shù)及其采樣值都變化很小，如果用連續(xù)的幀增加UV覆蓋顯然效果不好，為此需要隔一段時(shí)間增加UV覆蓋。UV覆蓋積分得到的是一段時(shí)間間隔完整數(shù)據(jù)幀疊加的數(shù)據(jù)，同樣對(duì)數(shù)據(jù)UU，VV及其采樣值操作得到的數(shù)據(jù)量成倍增加，這種操作能極大改善稀疏采樣圖像的質(zhì)量，提高信噪比。

圖4 UV覆蓋積分流程
Fig.4 Procedure of UV coverage overlap

以積分60 min的數(shù)據(jù)為例，積分后增加UV覆蓋，采樣數(shù)據(jù)成倍增加，對(duì)稀疏采樣提高圖像質(zhì)量及其信噪比起至關(guān)重要的作用。高頻陣單幀的UV覆蓋包含數(shù)據(jù)點(diǎn)1 770個(gè)，共軛之后有3 540個(gè)數(shù)據(jù)，得到一張2 560 × 2 560像素的圖像，覆蓋率僅僅0.054%。而對(duì)數(shù)據(jù)總共積分60 min間隔10 min疊加UV覆蓋后共有10 620個(gè)數(shù)據(jù)點(diǎn)，共軛之后數(shù)據(jù)點(diǎn)21 240個(gè)，覆蓋率提升為0.324%。選取60 ms采樣數(shù)據(jù)進(jìn)行平均，之后每隔10 min增加UV覆蓋，得到的UV覆蓋圖像如圖5、圖6。

圖5 3 ms單幀UV覆蓋
Fig.5 UV Coverage for a single frame (3ms)

圖6 60 min積分UV覆蓋
Fig.6 UV Coverage for 1 hour′s integration

2.2 積分代碼實(shí)現(xiàn)

為了實(shí)現(xiàn)UV平均積分、UV覆蓋積分以及UV平均和UV覆蓋混合積分，采用Python編寫(xiě)，同時(shí)為了簡(jiǎn)化程序，將3種模式的積分寫(xiě)入一個(gè)接口，提供模式的選擇，程序接口形式如下：

def clean_integration_R(self, sub_ARRAY, is_loop_mode, start_time, end_time, TASK_TYPE, time_average,

time_interval, BAND, CHANNEL, PLOT_ME, WRITE_FITS, P_ANGLE, DEBUG, outdir):

主要參數(shù)：sub_ARRAY為MUSER期數(shù)的選擇MUSER-I或者M(jìn)USER-II，start_time和end_time為開(kāi)始積分時(shí)間和結(jié)束積分時(shí)間，TASK_TYPE為積分模式的選擇包含UV平均積分、UV覆蓋積分、UV平均和UV覆蓋混合積分，time_average為UV平均積分的積分時(shí)間，time_interval為UV覆蓋積分時(shí)間間隔即間隔多長(zhǎng)時(shí)間疊加一次UV覆蓋。在系統(tǒng)中的運(yùn)行情況如圖7。代碼編寫(xiě)可以適用于數(shù)據(jù)采集的任何時(shí)間段積分模式。

圖7 MUSER系統(tǒng)運(yùn)行積分命令
Fig.7 Integration task in MUSER data processing System

3 實(shí)驗(yàn)分析

積分實(shí)例以混合積分模式為例，進(jìn)行10 min、20 min、30 min、40 min、50 min以及60 min的數(shù)據(jù)做積分處理，對(duì)60 ms的數(shù)據(jù)進(jìn)行UV平均積分，之后每隔10 min增加UV覆蓋。臟圖結(jié)果如圖8。

圖8從左到右從上到下分別為積分10，20，30,…,60 min之后的臟圖，可以很明顯地看出從第1張到最后一張?zhí)?yáng)輪廓略清晰。隨著積分時(shí)間的增加得到的臟圖質(zhì)量越來(lái)越好，圖像越來(lái)越清晰，同時(shí)可以清楚地看到太陽(yáng)的輪廓(圖中間類(lèi)似圓形)以及太陽(yáng)爆發(fā)的亮點(diǎn)位置，很明顯積分達(dá)到提升圖像質(zhì)量的效果。經(jīng)過(guò)潔化處理得到的潔圖如圖9。

圖9從左到右從上到下分別為積分10，20，30，…，60 min之后的潔圖，可以很明顯地看出，隨著積分時(shí)間的增加，得到的潔化圖像越來(lái)越清晰，同時(shí)可以清楚地看到太陽(yáng)的輪廓以及太陽(yáng)爆發(fā)的位置亮點(diǎn)，從一開(kāi)始多處亮點(diǎn)爆發(fā)到最后一張只有一處非常明顯的爆發(fā)點(diǎn)，信噪比有很大提升。

圖10為2016年7月5日04:04:38UT，射電日像儀高頻陣在4.187 8 GHz上用60 min積分時(shí)間觀測(cè)得到的太陽(yáng)射電圖像結(jié)果。與日本野邊山射電日像儀2016年7月5日04:00:02UT在17 GHz上觀測(cè)得到的太陽(yáng)圖像比較，高頻陣得到的太陽(yáng)像與野邊山射電日像儀觀測(cè)得到的太陽(yáng)像結(jié)果基本一致，可見(jiàn)積分方法的實(shí)現(xiàn)結(jié)果比較可靠。

圖8 混合積分臟圖Fig.8 Dirty images of mixed integration

圖9 混合積分潔圖
Fig.9 Clean images of mixed integration

圖10(a) MUSER-II 60 min積分的太陽(yáng)像；(b) 野邊山射電日像儀觀測(cè)的太陽(yáng)像

Fig.10(a) One huor′s integration solar images of MUSER-II; (b) Solar observation images of the Nobeyama Radioheliograph

4 結(jié) 論

通過(guò)對(duì)采樣數(shù)據(jù)積分處理，以每60 ms UV積分和間隔10 min增加UV覆蓋混合積分方法為例，可以看出成圖效果非常明顯，經(jīng)過(guò)一系列數(shù)據(jù)處理計(jì)算得到的臟圖可以清楚地看出太陽(yáng)邊緣輪廓及其爆發(fā)位置，信噪比有了很大提升。并且隨著積分時(shí)間的增加，得到的圖像質(zhì)量也有所提升。這一積分思路對(duì)于獲得質(zhì)量較優(yōu)的MUSER觀測(cè)結(jié)果有較好的作用。目前圖像處理結(jié)果得到很大提升，但在成圖前需要讀取UVFITS數(shù)據(jù)文件，當(dāng)前的UVFITS文件生成時(shí)間較慢，急需加快生成的效率。后續(xù)積分工作可以在積分時(shí)間上進(jìn)行探索，利用評(píng)價(jià)指標(biāo)尋找如何積分才能達(dá)到相對(duì)比較好的效果，給出UV積分時(shí)長(zhǎng)、UV覆蓋積分時(shí)間間隔以及總共60 min積分時(shí)間等參數(shù)的量化值。

[1] 向德琳. 射電天文觀測(cè)[J]. 北京: 科學(xué)出版社,1990.

[2] Gary D E, Keller C U. Solar and space weather radiophysics-current status and future developments[M]. Netherland: Kluwer Academic Publishers, 2004.

[3] Wang W, Yan Y, Liu D, et al. Calibration and data processing for a Chinese spectral radioheliograph in the decimeter wave range[J]. Publications of the Astronomical Society of Japan, 2013, 65(sp1): 2226-2237.

[4] 顏毅華, 張堅(jiān), 陳志軍, 等. 關(guān)于太陽(yáng)厘米-分米波段頻譜日像儀研究進(jìn)展[J]. 天文研究與技術(shù)——國(guó)家天文臺(tái)臺(tái)刊, 2006, 3(2): 91-98.

Yan Yihua, Zhang Jian, Chen Zhijun, et al. Progress on Chinese solar radioheliograph in cm-dm wavebandes [J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2006, 3(2): 91-98.

[5] 梅盈. MUSER海量數(shù)據(jù)預(yù)處理關(guān)鍵技術(shù)研究[D]. 昆明: 昆明理工大學(xué), 2015.

[6] Nakajima H, Nishio M, Enome S, et al. The Nobeyama radioheliograph[J]. Proceedings of the IEEE, 1994, 82(5): 705-713.

[7] Grechnev V V, Lesovoi S V, Smolkov G Y, et al. The Siberian Solar Radio Telescope: the current state of the instrument, observations, and data[J]. Solar Physics, 2003, 216(1-2): 239-272.

[8] Wang F, Mei Y, Deng H, et al. Distributed data-processing pipeline for Mingantu Ultrawide Spectral Radioheliograph[J]. Publications of the Astronomical Society of the Pacific, 2015, 127(950): 383-396.

[9] Wei S, Wang F, Deng H, et al. OpenCluster: a flexible distributed computing framework for astronomical data processing[J]. Publications of the Astronomical Society of the Pacific, 2016, 129(972): 1-14.

[10]Thompson A R, Moran J M, Swenson Jr G W. Interferometry and synthesis in radio astronomy[M]. Berlin: Springer, 2001: 648-648.

[11]羅洪禮. MUSER成像中的Grid技術(shù)研究與實(shí)現(xiàn)[D]. 昆明: 昆明理工大學(xué), 2016.

[12]Ségransan D. Observability and UV coverage[J]. New Astronomy Reviews, 2003, 51(8-9): 597-603.

[13]Mei Y, Liu D, Wang F, et al. The data format preprocessing system for Chinese spectral radioheliograph[J]. ICIC Express Letters, Part B: Applications, 2015, 6(5): 1467-1472.