2. 昆明理工大学, 云南 昆明 650051
2. Kunming University of Science and Technology, Kunming 650051, China
近十年来我国在射电天文领域得到了飞速发展，21 cm望远镜阵列(The 21 CentiMeter Array, 21CMA[1])、新疆的天籁阵列[2]、内蒙古的明安图太阳频谱日像仪(Mingantu Ultrawide Spectral Radioheliograph, MUSER) [3]以及中国甚长基线干涉测量(Very Long Baseline Interferometry, VLBI)网，这些射电望远镜的建成在我国天文科学研究、数据处理、经验积累、人才培养等方面都做出了重要贡献。平方千米阵(Square Kilometre Array, SKA)是人类正在建设的最大规模射电干涉望远镜，由13个成员国共同建设，中国是其中7个发起国之一。
观测数据的存储和交换是天文观测和研究的基本要求。FITS (Flexible Image Transport System, FITS)一直是天文数据存储交换的标准格式，针对射电观测数据，在FITS基础上发展出了UVFITS和FITSIDI等格式。近年来，测量集应用越来越广泛，逐渐成为射电天文数据处理分析的标准格式，被CASA、WSCLEAN[4]等射电天文数据处理软件广泛支持。测量集在国内射电天文领域应用相对较少，国内射电望远镜往往根据各自接收机的特点，自行定义相应的原始数据存储格式。如明安图太阳频谱日像仪采用祼二进制的方式保存观测文件，以大幅度降低存储空间，需要进行数据交换时，通过格式转换软件转换为UVFITS或FITSIDI。
射电天文模拟校准成像库是目前正在研制中的平方千米阵的算法参考库(Algorithm Reference Library, ARL)。为了实现射电天文模拟校准成像库与主流射电天文数据处理软件(如CASA)的数据对接，需要解决测量集格式输出问题。本文结合实际需求，系统讨论并实现了利用Python和Python-casacore生成测量集格式文件，并将Python代码集成到射电天文模拟校准成像库，为平方千米阵科学数据处理研究提供支撑，对其它射电望远镜数据处理工作也有较高的参考价值。
1 测量集文件测量集是一个遵从射电干涉测量方程(Radio Interferometer Measurement Equation, RIME[5])的文件格式标准，在AIPS ++ Note 191中被正式定义，用于规范校准前的射电天文观测数据的存储。测量集标准发布以后，CASA团队和欧洲VLBI网团队的多个天文软件开发小组进行了代码实现。由于CASA采用测量方程作为其基本校准方案，测量集很自然地成为CASA观测数据的存储标准。随着CASA成为阿塔卡玛大型毫米波天线阵和甚大阵(Very Large Array, VLA)的指定数据处理分析软件包[6]，测量集也成为它们数据分析中的缺省数据格式。但是，阿塔卡玛大型毫米波天线阵和甚大阵的原始数据存储格式分别为阿尔马科学数据模型(ALMA Science Data Model, ASDM)和科学数据模型(Science Data Model, SDM)，因此也都开发了相应接口，实现ASDM/SDM与测量集的转换。
测量集适用于目前射电天文学中的所有用例，包括单碟、少数天线构成的简单干涉仪以及成千上万个天线构成的大型射电干涉阵。测量集借鉴了关系型数据库的建模方法来降低数据的冗余，采用主表和子表、主键和外键，把干涉得到的可见度函数或单天线总功率测量值及其时间戳保存在主表(建主键)，反复使用的元数据保存在子表(建外键)，通过外键实现对元数据的引用。主表中有多个数据列以及关联到子表的键值，主表中必须有DATA列，用于存放干涉阵的可见度数据，或者FLOAT_DATA列，用于存放单天线功率值。
CASA遵循了测量集第2版本(MSv2) [7]，为了确保与CASA的数据兼容，本文在测量集第2版本基础上开展研究工作。
1.1 测量集的子表子表中存储了测量集的关键字，表 1列出了CASA采用的测量集第2版的所有子表，其中括号里的子表为可选子表。在实际应用中，必须生成非可选子表，即必须有的子表，子表内容可以为空。每个测量集文件一定要有一个主表(MAIN)，表中包含数据列和各子表的键。
| Table | Contents | Keys |
| ANTENNA | Antenna characteristics | ANTENNA_ID |
| DATA_DESCRIPTION | Data description | DATA_DESC_ID |
| (DOPPLER) | Doppler tracking | DOPPLER_ID, SOURCE_ID |
| FEED | Feed characteristics | FEED_ID, ANTENNA_ID, TIME, SPECTRAL_WINDOW_ID |
| FIELD | Field position | FIELD_ID |
| FLAG_CMD | Flag commands | TIME |
| (FREQ_OFFSET) | Frequency offset information | FEED_ID, ANTENNAn, FEED_ID, TIME, SPECTRAL_WINDOW_ID |
| HISTORY | History information | OBSERVATION_ID, TIME |
| OBSERVATION | Observer, Schedule, etc | OBSERVATION_ID |
| POINTING | Pointing information | ANTENNA_ID, TIME |
| POLARIZATION | Polarization setup | POLARIZATION_ID |
| PROCESSOR | Processor information | PROCESSOR_ID |
| (SOURCE) | Source information | SOURCE_ID, SPECTRAL_WINDOW_ID, TIME |
| SPECTRAL_WINDOW | Spectral window setups | SPECTRAL_WINDOW_ID |
| STATE | State information | STATE_ID |
| (SYSCAL) | System calibration characteristics | FEED_ID, ANTENNA_ID, TIME, SPECTRAL_WINDOW_ID |
| (WEATHER) | Weather info for each antenna | ANTENNA_ID, TIME |
显然，与FITS相比，测量集格式要复杂很多。在实际应用中，根据射电望远镜的不同观测数据需求，生成相应的子表并存入相关的数据，最终构成一个完整的测量集格式文件目录树结构，原则上必须生成所有非可选表。
与其它天文文件存储格式不同，测量集格式采用多级目录多文件保存的方法，各个表都以CASA表格存储。这意味着一个表格包含了多个文件，整个测量集格式文件也不是单个文件，而是一个由多级目录构成的目录树。一般来说，主表位于第1级目录，各个子表位于第2级目录。每一级目录均包含实际数据存放位置信息的table.info, table.f0, table.f1等文件。
1.2 测量集的主表结构每个测量集文件必须有一个主表(MAIN TABLE)，结构见表 2。受限于篇幅，表 2仅列出了主表的部分字段，完整的主表字段请参考测量集技术规范[7]。
| MAIN table: Data, coordinates and flags | ||||
| Name | Format | Units | Measure | Comments |
| Columns | ||||
| Keywords | ||||
| MS_VERSION | Float | MS format version | ||
| key | ||||
| TIME | Double | s | EPOCH | Integration midpoint |
| (TIME_EXTRA_PREC) | Double | s | extraTIME precision | |
| ANTENNA1 | Int | First antenna | ||
| …… | ||||
| DATA_DESC_ID | Int | Data desc. id. | ||
| PROCESSOR_ID | Int | Processor id. | ||
| FIELD_ID | Int | Field id. | ||
| Non-key attributes | ||||
| INTERVAL | Double | s | Sampling interval | |
| EXPOSURE | Double | s | The effective integration time | |
| …… | ||||
| OBSERVATION_ID | Int | Observation id. | ||
| STATE_ID | Int | State id. | ||
| UVW | Double(3) | m | UVW | UVW coordinates |
| Data | ||||
| (DATA) | Complex(Nc, Nf) | Complex visibility matrix (synthesis arrays) | ||
| (FLOAT_DATA) | Float(Nc, Nf) | Float data matrix (single dish) | ||
| Flag information | ||||
| FLAG | Bool(Nc, Nf*) | Cumulative data flags | ||
| 注：Nc为相关器数；Nf为频率通道数；Ncat为分类数；*表示不同的天线类型有不同的取值，对于本文讨论的射电干涉阵来说，都是相关参数的数量。 | ||||
由表 2可见，在数据存储方面，测量集格式与FITS格式基本类似，需要保存的数据类型也包括整型(int)、浮点(float)、双精度(double)、字符串(string)等。与FITS文件的头定义相比，测量集格式文件字段设计更为复杂，有3种类型的字段：
(1) 关键字(Keywords)：MS_VERSION，用来标识所保存的测量集格式文件遵从哪一个版本的规范。
(2) 键(Key)：相当于关系型数据库中的主键，用来和子表进行关联。如TIME给出了观测时刻，ANTENNA1，ANTENNA2指定两个相关天线。
(3) 非键属性(Non-key attributes)：根据实际需要定义的重要参数或属性。
其它表结构与主表类似，也各自规定了相应的保留字、键值和参数。生成测量集格式文件，必须首先明确各字段的内容、格式、单位等。
2 基于Python-casacore的测量集生成方法从AIPS++发展到CASA后，CASA软件的开发采用了多计算机语言混合编程的方法，其核心代码来自于C/C++，用户接口与调用部分基本上采用Python实现，C/C++部分构成了CASA核心库，即所谓的Casacore。Casacore是目前最完整的射电天文数据处理软件包，也是唯一实现了测量集文件的读写操作的软件包。Python-casacore https://github.com/casacore/casacore实现了Python对Casacore的调用。本文深入分析了Python-casacore的基本函数，研究其数据表的读写方法，完成了一个测量集格式输出对象，并集成到射电天文模拟校准成像库中，实现射电天文模拟校准成像库的测量集格式输出功能。
2.1 Python-casacore的表输出操作Python通过Python-casacore访问Casacore，表 3列出了Python-casacore的重要测量集表写入函数。
| No. | Function | Description |
| 1 | Table | Get the table object this column belongs to |
| 2 | putinfo | Put the table info. The table info is a dictionary containing the fields |
| 3 | putkeyword | Put the value of a table′s keyword |
| 4 | putkeywords | Put the value of multiple table keywords |
| 5 | putcol | Put an entire column or part of it |
| 6 | putcell | Put a value into one or more table cells |
| 7 | tableutil.makearrcoldesc | Create description of column holding arrays |
| 8 | tableutil.makescacoldesc | Create description of column holding scalars |
| 9 | tableutil.maketabdesc | Create table description |
| 10 | close | Flush and close the table which invalidates the table object |
调用表 3中相应的函数或者类，可以生成一个完整的数据子表，主要代码见表 4：
| No. | code | Description |
| 1 | col1 = tableutil.makearrcoldesc(‘UVW’, 0.0, 1, comment=‘Vector with uvw coordinates (in meters)’, keywords={‘QuantumUnits’: [‘m’, ‘m’, ‘m’], ‘MEASINFO’: {‘type’: ‘uvw’, ‘Ref’: ‘ITRF’}} |
Generate the column of col1 with fields of UVW, float type, meter unit, for storing a 1-dimension array |
| 2 | col2 = tableutil.makearrcoldesc(‘FLAG’, False, 2, comment=‘The data flags, array of bools with same shape as data’) |
Generate the column of col2 with the field of FLAG, boolean type, for storing a 2-dimension array |
| 3 | desc = tableutil.maketabdesc([col1, col2, col6, col7]) tb = table(“%s” % self.basename, desc, nrow=0, ack=False) |
Define a table desc with columns of col1, col2, col6, and col7. Generate an empty table tb using the definition of desc |
| 4 | for j in range(nBand): fg = numpy.zeros((nBL, self.nStokes, self.nchan), dtype=numpy.bool) tb.putcol(‘UVW’, uvwList, i, nBL) tb.putcol(‘FLAG’, fg.transpose(0, 2, 1), i, nBL) |
Write the data into the table tb |
(1) 生成矩阵列col1，定义‘UVW’字段，浮点类型，单位为m，该列存储一维数组。
(2) 生成矩阵列col2，定义‘FLAG’字段，布尔类型，该列存储二维数组。
(3) 生成表定义desc，该表有4列，分别为col1，col2，col6，col7。然后按照表定义desc创建空表tb。
(4) 把数据写入表tb，循环次数等于波段数量。
3 基于Python-casacore的测量集生成及其与RASCIL的集成应用本文利用Python-casacore，针对测量集第2版开发实现了测量集文件输出模块，封装成WriteMS类，结构见图 1。
|
| 图 1 WriteMs类图 Fig. 1 The class of WriteMs |
目前，开发的程序已经集成到平方千米阵的算法参考库RASCIL (https://gitlab.com/timcornwell/rascil/tree/master/rascil/processing_components/visibility)，同时，在射电天文模拟校准成像库中给出了一个利用WriteMS输出测量集格式的例子，供使用者参考，见test_export_ms_rascil.py。图 2是该应用实例的流程图，包括了模拟观测到测量集输出的完整过程，主要步骤解释如下：
|
| 图 2 测量集文件生成应用实例流程图 Fig. 2 The flow chart of exporting MS format file |
(1) Input parameters。输入参数主要包括观测位置、观测时间、观测频率、频带宽度、相位中心、输出文件的名称以及干涉阵列天线配置。
(2) Input the image of M31。实例中使用了阿塔卡玛大型毫米波天线阵和其它模拟程序通常采用的M31星系图作为模拟观测图像，见图 3(a)。
|
| 图 3 M31的原始图像与模拟观测图像的比较。(a)原始图像；(b)模拟观测脏图 Fig. 3 The original M31 image and the observed M31 image. (a) The original image; (b)The simulated image |
(3) Create block of visibilities。根据步骤(1)所输入的参数和天线配置生成初始可见度数组，预测宽视场成像参数，得到合适的栅格尺寸大小。
(4) Sample M31。对M31图像进行模拟采样，得到观测可见度数组。
(5) Export visibilities to MS files。把重采样得到的可见度函数按测量集格式输出。
为了检验本文方法的正确性，使用CASA软件对所生成的测量集文件成图。即把图 2流程所生成的测量集格式的可见度数据导入CASA，利用CASA成图(图 3(b))，可以看出所成的脏图与原始图像(图 3(a))轮廓基本一致。因为CASA软件在读取测量集格式数据的过程中会进行较多的校验操作，因此，如果所生成的测量集格式数据通过了CASA的处理，得到了正确的结果，说明所生成的测量集各个子表和字段是合理的(虽然部分字段的取值可能与真实情况有区别，但这不影响成像的处理)。
4 结论虽然测量集文件规范制定较早，但在我国的射电领域应用较少。一方面是因为测量集格式文件占用较多的空间，另一方面是因为生成测量集格式文件一直依赖于Casacore这一底层软件包，开发比较困难。本文结合平方千米阵工程建设的需要，研究了测量集格式文件的定义、结构、字段以及利用Python-casacore写入数据的方法，并基于Python-casacore开发实现了测量集格式输出软件包，集成到射电天文模拟校准成像库。整体来看，本文的工作在平方千米阵算法参考库的研制过程中起到了重要的支撑作用，为后续数据模拟与文件存储提供了保障，也对其它射电天文数据的测量集格式文件生成有较好的参考作用。
致谢: 感谢国家天文台-阿里云天文大数据联合研究中心对本项工作的支持。
| [1] | WU X P. Probing the epoch of reionization with 21CMA: status and prospects[C]//Proceedings of Science. 2009. |
| [2] | 陈学雷. 暗能量的射电探测——天籁计划简介[J]. 中国科学:物理学力学天文学, 2011, 41(12): 1358–1366 |
| [3] | YAN Y H, ZHANG J, HUANG G L. On the Chinese spectral radioheliograph (CSRH) project in cm-and dm-wave range[C]//2004 Asia-Pacific Radio Science Conference. 2004: 391-392. |
| [4] | OFFRINGA A, MCKINLEY B, HURLEY-WALKER N, et al. WSCLEAN:an implementation of a fast, generic wide-field imager for radio astronomy[J]. Monthly Notices of the Royal Astronomical Society, 2014, 444(1): 606–619. |
| [5] | HAMAKER J, BREGMAN J, SAULT R. Understanding radio polarimetry. I. mathematical foundations[J]. Astronomy and Astrophysics Supplement Series, 1996, 117(1): 137–147. |
| [6] | PETRY D. Analysing ALMA data with CASA[C]//ASP Conference Series. 2012. |
| [7] | KEMBALL A, WIERINGA M H. MeasurementSet definition version 2.0[EB/OL]. 2000[2019-08-08]. https://casa.nrao.edu/Memos/229.html. |