MUSER可见度数据积分方法与实现
赖铖1, 梅盈1,2, 邓辉1, 王锋1,2, 戴伟1     
1. 昆明理工大学云南省计算机技术应用重点实验室, 云南 昆明 650500;
2. 中国科学院云南天文台, 云南 昆明 650011
摘要: 射电观测是研究太阳活动的重要探测手段。我国明安图射电频谱日像仪(MingantU SpEctral Radioheliograph,MUSER)主要用于研究太阳爆发活动初始能量释放区的物理过程,其观测将在太阳射电成像开辟一个新的窗口。成像处理是数据处理的重要组成部分,如何提高成像质量是当前数据处理的研究重点。首先介绍了射电干涉成像的基本理论,随后分析了对观测得到的可见度数据积分的必要性,细致讨论了短时段可见度数据叠加求平均和长时段UV覆盖叠加两种积分方法,并给出了完整的实现。通过实现代码与实验验证,两种积分均可以有效提高信噪比,图像质量明显提高。
关键词: MUSER     UV平均积分     UV覆盖积分     综合孔径成像    
Integral Method and Implementation of MUSER Visibility Data
Lai Cheng1, Mei Ying1,2, Deng Hui1, Wang Feng1,2, Dai Wei1     
1. Key Laboratory of Applications of Computer Technology of the Yunnan Province, Kunming University of Science and Technology, Kunming 650500, China;
2. Yunnan Observatories, Chinese Academy of Science, Kunming 650011, China
Abstract: Radio observation is the most important detection means to research the solar activity. The MingantU SpEctral Radioheliograph (MUSER) is mainly used to research the physical process in the initial energy-releasing area during the solar activity and its observation will open a new window in solar radio imaging. Imaging processing is an important part of MUSER data processing and improving the image quality is the most important research of data processing. Firstly, we introduce the basic theory of radio interference imaging. Secondly, we analyze the necessity of visibility data integration obtained from the MUSER observations. Then, we discuss in detail the two kinds of integration method which are short-term visibility data superimposed averaging and long-term UV overlay stacking and give out a complete realization. Through the method of verification and the validity of completion, the two kinds of integration method can effectively improve the signal-to-noise ratio and the quality of images.
Key words: MUSER     UV integration     UV coverage integration     Aperture synthesis imaging    

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

MUSER天线阵包含低频阵(MUSER-Ⅰ)和高频阵(MUSER-Ⅱ)。低频阵共有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覆盖积分两种方法提高太阳图像的信噪比。

http://secchirh.obspm.fr/index.php

1 成像原理 1.1 综合孔径成像原理

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

$ V\left( {u, v} \right) = \sum \sum I\left( {l, m} \right){e^{{\rm{i}}2{\rm{ \mathsf{ π} }}(ul + vm)}}. $ (1)

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

$ V\left( u, v \right)=\iint{I\left( l, m \right){{e}^{{\rm{i2 }\rm{ \mathsf{ π}}}(ul+vm)}}}{\rm{d}}l{\rm{d}}m\rm{, } $ (2)

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

对于任意一个亮度分布I(l, m),总有一个复可见度函数V(u, v)与之相对应。复可见度函数是连续的复函数,而天线阵只能构成有限条基线,即只能获得复可见度函数的有限个采样点。采样函数$ S\left( {u, v} \right) = \sum\limits_k {\delta (u - {u_k}, v - {v_k})} $,其中,k为取遍所有(u, v)的点。采样函数的输出即为天线得到的真正数据,对这些数据进行傅里叶逆变换得到脏图ID(l, m)[11]

$ {{I}_{\rm{D}}}\left( l, m \right)=\iint{\mathit{S}}{\left( u, v \right)}V\left( u, v \right){{e}^{-{\rm{i2 }} {\rm{ \mathsf{ π}}}(ul+vm)}}{\rm{d}}u{\rm{d}}v. $ (3)
1.2 MUSER积分的基本原理

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

图 1 MUSER原始脏图(a)和洁化图(b) Figure 1 Dirty and clean image without integration (MUSER-Ⅱ, 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积分表示如下:

$ \begin{align} & {\rm{采样函数:}}\ \ {{S}_{{{t}_{I}}}}\left( u, v \right)=\sum\limits_{i=1}^{\mathit{N}}{{{S}_{{{t}_{i}}}}({{u}_{i}}, {{v}_{i}})}\ \ \ (N=\frac{{{t}_{I}}}{{{t}_{c}}}). \\ & {\rm{采样输出:}}\ \ {{V}_{{{t}_{I}}}}\left( u, v \right)=\sum\limits_{i=1}^{\mathit{N}}{{{V}_{{{t}_{i}}}}({{u}_{i}}, {{v}_{i}})}\ \ \ (N=\frac{{{t}_{I}}}{{{t}_{c}}})\rm{ }. \\ & {\rm{UV平均脏图:}}\ \ {{I}_{\rm{D}}}\left( l, m \right)=\iint{\left[\frac{{{S}_{{{t}_{I}}}}\left( u, v \right)}{N} \right]}\left[\frac{{{V}_{{{t}_{I}}}}\left( u, v \right)}{N} \right]{{e}^{{-\rm{i2 } \rm{ \mathsf{ π}}}(ul+vm)}}{\rm{d}}u{\rm{d}}v\ \ (N={{\frac{{{t}_{I}}}{{{t}_{c}}}}}). \\ \end{align} $ (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覆盖后的采样函数和采样输出表示为

$ \begin{align} & {\rm{采样函数:}}\ \ {{S}_{{{t}_{J}}}}\left( u, v \right)=\sum\limits_{j=1}^{N}{{{S}_{{{t}_{j}}}}({{u}_{j}}, {{v}_{j}})}\ \ \ \ (N=\frac{{{t}_{a}}}{{{t}_{{\rm{int}}erval}}}). \\ & {\rm{采样输出:}}\ \ {{V}_{{{t}_{J}}}}\left( u, v \right)=\sum\limits_{j=1}^{N}{{{V}_{{{t}_{j}}}}({{u}_{j}}, {{v}_{j}})\ }\ \ \ (N=\frac{{{t}_{a}}}{{{t}_{{\rm{int}}erval}}}). \\ & {\rm{UV覆盖增加脏图:}}\ \ {{I}_{\rm{D}}}\left( l, m \right)=\iint{{}}{{S}_{{{t}_{J}}}}\left( u, v \right){{V}_{{{t}_{J}}}}\left( u, v \right){{e}^{{-\rm{i2 } \rm{ \mathsf{ π}}}(ul+vm)}}{\rm{d}}u{\rm{d}}v. \\ \end{align} $ (5)

由(5)式可知,在间隔tinterval增加UV覆盖,总积分时间ta可以增加$ \frac{{{t_a}}}{{{t_{{\rm{int}}erval}}}}$倍的UV数据点,假设有n面天线构成的天线阵列,可以得到n(n-1)/2个数据采样点,通过增加UV覆盖,采样数据点总共有$ \left( {\frac{{{t_a}}}{{{t_{{\rm{int}}erval}}}}} \right)\left( {\frac{{\mathit{n}\left( {\mathit{n}{\rm{ + 1}}} \right)}}{2}} \right)$个,这对稀疏采样的成图质量有很大改善。

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

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

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

$ \begin{align} & {\rm{采样函数:}}\ \ {{S}_{\rm{mix}}}\left( u, v \right)=\sum\limits_{j=1}^{N}{\frac{\sum\limits_{i=1}^{M}{{{S}_{{{t}_{ij}}}}({{u}_{ij}}, {{v}_{ij}})}}{M}}\ \ \ \ \ (M=\frac{{{t}_{I}}}{{{t}_{c}}}, N=\frac{{{t}_{a}}}{{{t}_{{\rm{int}}erval}}}). \\ & {\rm{采样输出:}}\ \ {{V}_{\rm{mix}}}\left( u, v \right)=\sum\limits_{j=1}^{N}{\frac{\sum\limits_{i=1}^{M}{{{V}_{{{t}_{ij}}}}({{u}_{ij}}, {{v}_{ij}})}}{M}\ \ \ \ \ }(M=\frac{{{t}_{I}}}{{{t}_{c}}}, N=\frac{{{t}_{a}}}{{{t}_{{\rm{int}}erval}}}). \\ & {\rm{脏图:}}\ \ {{I}_{D}}\left( l, m \right)=\iint{{{S}_{\rm{mix}}}\left( u, v \right)}{{V}_{\rm{mix}}}\left( u, v \right){{e}^{{-\rm{i2 } \rm{ \mathsf{ π} }}(ul+vm)}}{\rm{d}}u{\rm{d}}v. \\ \end{align} $ (6)

在总时间ta内UV覆盖增加了$ {\frac{{{t_a}}}{{{t_{{\rm{int}}erval}}}}}$倍,每组增加的UV覆盖是在时间段tI对UV平均后再叠加的覆盖。

2 积分过程实现

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

图 2 MUSER积分成图过程 Figure 2 The procedure of MUSER integration
2.1 积分操作流程

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

图 3 可见度数据平均的计算 Figure 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覆盖积分流程 Figure 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覆盖 Figure 5 UV Coverage for a single frame (3ms)
图 6 60 min积分UV覆盖 Figure 6 UV Coverage for 1 hour′s integration
2.2 积分代码实现

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

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

图 7 MUSER系统运行积分命令 Figure 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 混合积分脏图 Figure 8 Dirty images of mixed integration

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

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

图 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上观测得到的太阳图像比较,高频阵得到的太阳像与野边山射电日像仪观测得到的太阳像结果基本一致,可见积分方法的实现结果比较可靠。

图 10 (a) MUSER-Ⅱ 60 min积分的太阳像;(b)野边山射电日像仪观测的太阳像 Figure 10 (a) One huor′s integration solar images of MUSER-Ⅱ; (b) Solar observation images of the Nobeyama Radioheliograph
4 结论

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

参考文献
[1] 向德琳. 射电天文观测[M]. 北京: 科学出版社, 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. http://cdmd.cnki.com.cn/Article/CDMD-10674-1015636038.htm
[6] Nakajima H, Nishio M, Enome S, et al. The Nobeyama radioheliograph[J]. Proceedings of the IEEE, 1994, 82(5): 705–713. DOI: 10.1109/5.284737
[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. DOI: 10.1086/680852
[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. http://cdmd.cnki.com.cn/Article/CDMD-10674-1016200370.htm
[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.
由中国科学院国家天文台主办。
0

文章信息

赖铖, 梅盈, 邓辉, 王锋, 戴伟
Lai Cheng, Mei Ying, Deng Hui, Wang Feng, Dai Wei
MUSER可见度数据积分方法与实现
Integral Method and Implementation of MUSER Visibility Data
天文研究与技术, 2018, 15(1): 78-86.
Astronomical Research and Technology, 2018, 15(1): 78-86.
收稿日期: 2017-05-15
修订日期: 2017-07-03

工作空间