MUSER可见度数据积分方法与实现*

赖 铖，梅 盈,2，邓 辉，王 锋,2，戴 伟

(1. 昆明理工大学云南省计算机技术应用重点实验室，云南 昆明 650500；2. 中国科学院云南天文台，云南 昆明 650011)

射电天文是天文学的一个分支，它使用射电望远镜系统在无线电波段研究来自深空(包括各类天体)的射电波[1]。为推动射电观测技术的发展，填补国际上对太阳耀斑能量初始释放区分米波段高分辨射电成像观测的空白[2-3]，我国已建成在0.4～15 GHz范围内的厘米-分米波日像仪——明安图射电频谱日像仪。它可以在超宽频带下同时以高时间、高空间和高频率分辨率观测太阳的动力学性质，探索太阳剧烈活动的起源[4]。

MUSER天线阵包含低频阵(MUSER-I)和高频阵(MUSER-II)。低频阵共有40面4.5 m天线，观测频率为400 MHz～2 GHz。低频阵的观测分为4个波段，两种极化方式(左旋、右旋)，其中每个波段包含16个通道。高频阵共60面2 m天线，观测频率范围为2 GHz～15 GHz，在33个波段和两种极化方式上观测，每个波段包含16个通道[5]。高频阵和低频阵每3 ms产生一帧观测数据(一个波段、一个极化方向上16个通道的数据)。如何将原始观测的可见度数据经过一系列校准和计算得到用于科学研究的太阳图像，是数据处理的关键。

目前，国外用于太阳射电成像观测的设备主要有日本野边山的太阳射电日像仪，积分时间为25 ms[6]；法国南茜米波太阳射电日像仪，积分时间1 s①；西伯利亚太阳射电日像仪，积分时间0.336 s[7]。这些射电日像仪积分处理过程一般直接在接收机完成，因此时间分辨率较低。我国位于内蒙古明安图观测站的明安图射电频谱日像仪已经建成，满足实时和事后科学数据处理的软件已开发完成[8-9]。为进一步提高成图质量，本文通过对成图基本原理的分析，进一步提出基于UV平均积分和UV覆盖积分两种方法提高太阳图像的信噪比。

1 成像原理

1.1 综合孔径成像原理

相关干涉仪是将间隔一定距离的两面天线接收的天体辐射分别传送到相关器中进行相关处理。由干涉成像原理可知，干涉仪的输出响应和天空亮度存在傅里叶变换关系，假设天空亮度分布为I(l,m)，干涉仪输出为V(u,v)，则有如下关系[10]：

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

(1)

将离散求和转换成连续积分可以得出如下公式[10]：

(2)

其中，u,v为投影平面也称作空间频率；l,m为平时的天空面。 (2)式表明干涉仪的输出是亮度分布傅里叶变换的一个(u,v)分量，随着u,v的不同可以得到不同空间频率分量的值，所有的空间频率分量总称复可见度函数。对复可见度函数进行傅里叶逆变换得到天空亮度分布，即图像，这被称为综合孔径成像原理[10]。

(3)

1.2 MUSER积分的基本原理

当前数据处理系统流水线已完成从原始观测数据格式到FITS文件格式的转换，观测数据校准到脏图以及洁化的成图过程[8]。以高频阵在2016年7月5日4时4分38秒的观测数据为例，生成脏图(图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))

可以看出，对单帧采样数据成像后的太阳图像信噪比很低，不能清楚地看出太阳爆发的位置及其轮廓。脏图洁化后得到的洁图，同样也不能清楚地分析出爆发位置。为提高日像仪重建的太阳图像的信噪比，得到质量较好的成像结果，需要对采样数据进行积分处理，改善图像质量，提高信噪比[12]。本文根据科学研究的需要，以及积分时段要求，从两方面分析，研究了UV平均积分(可见度数据求和再求平均)、UV覆盖积分(利用地球自转增加UV覆盖)两种方法。

1.2.1 基于UV平均的积分

明安图射电频谱日像仪选择观测模式为快照模式，积分时间为3 ms进行观测。如此高的时间分辨率情况下，为了提高空间分辨率，达到改善图像质量的效果，采取对3 ms观测得到的UV及其输出数据平均的方法实现。

以高频阵为例，其中单帧积分采样时间tc，假设对UV积分时间为tI，单帧采样函数S(u,v)，采样输出为V(u,v)，则UV积分表示如下：

(4)

(4)式做法即将采样的复可见函数进行平均，同时采样函数也进行平均，用平均采样函数和复可见度函数的方法实现UV平均积分。

1.2.2 基于UV覆盖增加的积分

UV覆盖对成图起至关重要的作用，由于采样函数S(u,v)的傅里叶变换即得到脏束，也称综合束，而脏束又直接影响成图，因此可通过增加UV覆盖提高图像质量。增加UV 覆盖有3种方式，(1)增加阵列望远镜个数；(2)增加子孔径有效半径；(3)利用地球自转[12]。增加阵列望远镜个数对成本有直接需求，在成本固定的情况下不考虑，增加子孔径有效半径可能导致采样率不足。在现有成本下通过地球自转增加UV覆盖是一种非常直接可行的方法。

随着地球自转UV发生改变，而短时间内UV的变化很小，因此假设每间隔tinterval增加UV覆盖，总积分时间ta，则增加UV覆盖后的采样函数和采样输出表示为

(5)

上述两种方式均可以提高成图的质量，将两种方法混合到一起同样可以提升图像质量，由此提出先进行UV平均积分再增加UV覆盖积分的方法。

1.2.3 基于UV平均和增加覆盖的混合积分

假设UV积分时间为tI，可以知道在该段时间内的采样函数为StI(ui,vi)，采样输出VtI(ui,vi)，每间隔tinterval增加UV覆盖，总积分时间ta，混合积分采样函数及其采样输出可以表示为

(6)

2 积分过程实现

明安图射电频谱日像仪的成图过程：目标源放射出天体辐射，然后干涉仪对其进行采样，采样得到的原始数据经过预处理相关校准[13]、对应极化和通道，之后进行积分处理，再网格化傅里叶逆变换得到脏图，最后经过洁化算法处理得到洁图，过程如图2。

图2 MUSER积分成图过程
Fig.2 The procedure of MUSER integration

2.1 积分操作流程

明安图射电频谱日像仪的每3 ms接收一帧观测数据，对这些数据进行积分操作得到积分数据，积分数据是观测得到并且经过预处理相关校准、对应极化和通道的可见度数据。UV平均积分具体操作流程如图3。

图3 可见度数据平均的计算
Fig.3 Calculation of the average of visibilities

高频阵中每一帧完整的数据包含33小帧。UV积分对连续的每一帧数据中UU，VV以及采样得到的实部、虚部数据进行累加再求平均，得到的UV平均积分数据量和单帧大小数量一样多，例如在高频阵中采样函数UV点为3 540个采样得到的实部虚部3 540个，则经过UV平均积分得到的数据同样是3 540个UV点及其采样值。

实现UV覆盖积分具体操作流程如图4。对高频阵中数据进行UV覆盖积分同样用完整的数据帧积分，不同的是由于短时间内采样函数及其采样值都变化很小，如果用连续的帧增加UV覆盖显然效果不好，为此需要隔一段时间增加UV覆盖。UV覆盖积分得到的是一段时间间隔完整数据帧叠加的数据，同样对数据UU，VV及其采样值操作得到的数据量成倍增加，这种操作能极大改善稀疏采样图像的质量，提高信噪比。

图4 UV覆盖积分流程
Fig.4 Procedure of UV coverage overlap

以积分60 min的数据为例，积分后增加UV覆盖，采样数据成倍增加，对稀疏采样提高图像质量及其信噪比起至关重要的作用。高频阵单帧的UV覆盖包含数据点1 770个，共轭之后有3 540个数据，得到一张2 560 × 2 560像素的图像，覆盖率仅仅0.054%。而对数据总共积分60 min间隔10 min叠加UV覆盖后共有10 620个数据点，共轭之后数据点21 240个，覆盖率提升为0.324%。选取60 ms采样数据进行平均，之后每隔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 积分代码实现

为了实现UV平均积分、UV覆盖积分以及UV平均和UV覆盖混合积分，采用Python编写，同时为了简化程序，将3种模式的积分写入一个接口，提供模式的选择，程序接口形式如下：

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):

主要参数：sub_ARRAY为MUSER期数的选择MUSER-I或者MUSER-II，start_time和end_time为开始积分时间和结束积分时间，TASK_TYPE为积分模式的选择包含UV平均积分、UV覆盖积分、UV平均和UV覆盖混合积分，time_average为UV平均积分的积分时间，time_interval为UV覆盖积分时间间隔即间隔多长时间叠加一次UV覆盖。在系统中的运行情况如图7。代码编写可以适用于数据采集的任何时间段积分模式。

图7 MUSER系统运行积分命令
Fig.7 Integration task in MUSER data processing System

3 实验分析

积分实例以混合积分模式为例，进行10 min、20 min、30 min、40 min、50 min以及60 min的数据做积分处理，对60 ms的数据进行UV平均积分，之后每隔10 min增加UV覆盖。脏图结果如图8。

图8从左到右从上到下分别为积分10，20，30,…,60 min之后的脏图，可以很明显地看出从第1张到最后一张太阳轮廓略清晰。随着积分时间的增加得到的脏图质量越来越好，图像越来越清晰，同时可以清楚地看到太阳的轮廓(图中间类似圆形)以及太阳爆发的亮点位置，很明显积分达到提升图像质量的效果。经过洁化处理得到的洁图如图9。

图9从左到右从上到下分别为积分10，20，30，…，60 min之后的洁图，可以很明显地看出，随着积分时间的增加，得到的洁化图像越来越清晰，同时可以清楚地看到太阳的轮廓以及太阳爆发的位置亮点，从一开始多处亮点爆发到最后一张只有一处非常明显的爆发点，信噪比有很大提升。

图10为2016年7月5日04:04:38UT，射电日像仪高频阵在4.187 8 GHz上用60 min积分时间观测得到的太阳射电图像结果。与日本野边山射电日像仪2016年7月5日04:00:02UT在17 GHz上观测得到的太阳图像比较，高频阵得到的太阳像与野边山射电日像仪观测得到的太阳像结果基本一致，可见积分方法的实现结果比较可靠。

图8 混合积分脏图Fig.8 Dirty images of mixed integration

图9 混合积分洁图
Fig.9 Clean images of mixed integration

图10(a) MUSER-II 60 min积分的太阳像；(b) 野边山射电日像仪观测的太阳像

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

4 结 论

通过对采样数据积分处理，以每60 ms UV积分和间隔10 min增加UV覆盖混合积分方法为例，可以看出成图效果非常明显，经过一系列数据处理计算得到的脏图可以清楚地看出太阳边缘轮廓及其爆发位置，信噪比有了很大提升。并且随着积分时间的增加，得到的图像质量也有所提升。这一积分思路对于获得质量较优的MUSER观测结果有较好的作用。目前图像处理结果得到很大提升，但在成图前需要读取UVFITS数据文件，当前的UVFITS文件生成时间较慢，急需加快生成的效率。后续积分工作可以在积分时间上进行探索，利用评价指标寻找如何积分才能达到相对比较好的效果，给出UV积分时长、UV覆盖积分时间间隔以及总共60 min积分时间等参数的量化值。

[1] 向德琳. 射电天文观测[J]. 北京: 科学出版社,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] 颜毅华, 张坚, 陈志军, 等. 关于太阳厘米-分米波段频谱日像仪研究进展[J]. 天文研究与技术——国家天文台台刊, 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海量数据预处理关键技术研究[D]. 昆明: 昆明理工大学, 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技术研究与实现[D]. 昆明: 昆明理工大学, 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.