德令哈13.7 m望远镜谱线OTF观测系统∗

孙继先逯登荣杨戟苏扬张少博周鑫林镇辉

(1中国科学院紫金山天文台南京210034)

(2中国科学院射电天文重点实验室南京210034)

1 引言

OTF(On-The-Fly)观测模式是一种应用于单天线成图观测的技术[1],观测时天线以恒定的速率对指定的区域进行扫描,多个采样点(简称ON点)测量数据共用1个参考点测量数据和黑体测量数据进行校准,后端数据被频繁读取以避免波束污染.OTF观测方法广泛应用到射电望远镜观测当中,如NRAO(National Radio Astronomy Observatory,美国国立射电天文台)12 m望远镜[1]、Nobeyama 45 m和ASTE 10 m望远镜[2].国内13.7 m望远镜于2011年首次使用该方法,2015年天马望远镜运用了OTF观测方法[3].多数望远镜OTF系统应用在单波束接收机,用于多波束接收机的较少,多波束接收机下实现OTF观测更为复杂,需要测试多个波束之间的相对位置以及各个波束的相对效率.接收机后端又有连续谱和谱线观测模式之分,连续谱观测通常只有1个通道,而谱线观测要观测频谱,有很多通道,如13.7 m望远镜后端频谱仪有16384个通道,后端数据采集量比连续谱观测大16384倍,从而对数据采集和存储有更高的要求.

13.7 m望远镜是典型的地平式卡塞格林结构天线,天线伺服系统有独立的控制系统,接收调度机的指令对天线进行控制.超导成像频谱仪[4−5]采用边带分离的超导混频器,有3×3个波束(beam),有边带分离SIS(Superconductor-Insulator-Superconductor)混频器、数字偏置源、独立中频、快速傅里叶变换数字频谱仪(FFTS)等主要部分,可以同时接收18路信号.对安装在地平式望远镜上的多波束接收机来说,中心波束跟踪天空的一个位置,由于像场旋转问题,其他波束对应的天空位置会随着时间改变,所以其他波束无法长时间跟踪固定位置进行积分.为了消除像场旋转,必须配备专门的像场消转装置,由于接收机比较重,像场消转装置会增加系统复杂程度.因为OTF观测频繁读取后端数据和天线位置,通过后期数据处理,可以计算出各个波束指向的天空位置,经过网格化(Gridding)处理就能得到规则的成图数据,所以不用像场消转装置就可以实现多波束的成图观测.

本文介绍13.7 m望远镜的OTF观测系统,以及具体观测过程,包括系统软硬件组成、观测过程和参数优化、数据预处理相关参数测试及处理过程等.

2 OTF观测系统的软件和硬件

2.1 系统结构

13.7 m望远镜控制系统是分布式的,OTF观测涉及到很多分系统,有调度计算机、天线控制单元(ACU)、副面控制单元(SRCU)、前端控制单元(FECU)、后端控制单元(BECU)、两台数据采集单元(FFTS)和数据预处理系统(DPS).OTF控制软件协调各个单元完成观测过程,各个单元之间通讯采用网络通讯.每个OTF观测扫描结束时在记录文件中追加一条记录,每小时往DPS上传1次记录.DPS上设置定时任务每小时执行1次,根据记录对数据进行实时预处理,并将预处理后的CLASS格式和FITS格式数据上传到毫米波射电天文数据库,供观测者下载.

2.2 控制软件

OTF控制软件在调度计算机上运行,运行环境采用Linux系统,程序用C语言编写.控制软件包括监控界面程序(Monitor)、天线控制程序和OTF观测程序3个进程,3个进程分别在3个终端独立运行,进程之间采用共享内存的方式实现相互通信.Monitor进程显示时间、天线、接收机和观测源等信息.天线控制程序根据观测程序指令计算天线的位置及运转速度,并将指令发送给ACU,每50 ms读取一次天线位置、状态等信息.OTF观测模式和其他观测模式(如单点观测、指向观测等)共用Monitor和天线控制程序,切换观测模式只需改变观测程序即可.OTF观测程序从星表文件中读取观测源的坐标和OTF扫描参数,观测过程中根据需要给FECU和BECU发指令控制接收机状态.每扫描一行读取一次接收机状态和天气信息,如果接收机不正常或有云,该行的数据剔除标志标记为1.观测黑体时给SRCU发指令,让副面进行自适应控制[6−8].天线控制程序根据观测进程给ACU发指令让天线跟踪到需要的位置,当天线跟踪到位后观测程序给FFTS发指令开始采集数据.

2.3 数据采集

数据采集系统为基于FPGA(Field Programmable Gate Array)的宽带、高分辨FFTS.每台频谱仪安装有9块采集板卡,接受9路独立的IF(Intermediate Frequency)输入.每路独立进行8 bit、2 GHz采样,瞬时带宽为1 GHz.采样后的数据流用FPGA进行实时快速傅里叶变换(FFT)频谱运算,信号带宽为1 GHz,频点数为16384,频率分辨率为61 kHz.

由两台FFTS分别采集上下边带的数据,FFTS操作系统为Linux,程序采用C语言编写.9个子线程独自从FPGA读取数据,主线程协调采集,数据采集校准后保存到数据预处理系统,每个波束扫描一遍保存一个文件.数据预处理系统挂载了40 TB的磁盘阵列,通过nfs共享方式将磁盘共享给FFTS.

FFTS采集的数据校准成温度后以二进制格式保存,每行扫描结束时保存一次.每行的数据为一块,前面是观测源的信息(如坐标、谱线名称、视向速度、观测频率)、波束号、FFTS采样每个点时的计算机时间、中心通道数、日期、一行采样点数、此行数据是否删除标志以及每个采样点对应的地方恒星时(local sidereal time,LST)、儒略日、方位角(AZ)、仰角(EL)等,后面部分是该行所有观测点的谱线数据.

3 OTF观测过程和参数优化

3.1 制作观测星表

观测之前先要制作观测星表,星表中包含源坐标和OTF扫描设置所需的参数.利用“德令哈OTF观测星表生成器”(DOTG,界面如图1)制作观测星表.DOTG用TCL/TK语言编写,具有友好的界面,在相应输入框中输入源名、视向速度、源中心点坐标、参考位置坐标.由于9个波束的波束间隔为3′,加上像场旋转,参考位置要求中心点周围±5′范围内没有发射谱线.该软件具有参考位置查找功能,用哥伦比亚1.2 m望远镜CO(2-1)分子谱线巡天的数据[9](数据下载网址https://skyview.gsfc.nasa.gov/current/cgi/moreinfo.pl?survey=CO),将积分强度小于3 K·km·s−1的数据作为参考星表,软件从参考星表查找在离源2°范围内的所有坐标,如果2°范围内没有找到,自动将查找范围增加1°再查找,最大范围为6°.自动选择离观测中心最近的积分强度小于3 K·km·s−1的位置,查找结果绘到图上,观测者可以在图上选择合适的位置作为参考位置,软件还绘出源和参考位置在方位角和仰角上的差值随LST变化的图,方便观测者选择仰角上差值小的时间段进行观测,这样能保证有比较平坦的基线.

扫描方式有赤道方向、银道方向和与赤道成某一夹角3种方式可选,其他OTF扫描参数都从界面输入.界面上默认的参数为经过测试比较理想的参数,观测者只需要修改扫描范围即可,软件会计算该参数下观测需要的时间和网格化为30′′间隔后能达到的噪声水平.

图1 德令哈OTF观测星表生成器软件界面Fig.1 Software interface of the Delingha OTF observational catalog generator

3.2 观测流程

观测过程中相关计算公式如下:

式中θRA为起始位置赤经(RA)偏移量,θx为扫描区域宽度,θramp为加速区域大小,c是一个函数,如果扫描行数i是偶数,结果为−1,i是奇数,则结果为1.

式中θDEC为起始位置赤纬(DEC)偏移量,θy为扫描区域高度,θRow为相邻两个扫描行的间隔.

式中p为星位角,ϕ为地理纬度,δ为赤纬,t为LST与RA的差值.

式中θAZ为起始位置换算到AZ上的偏移量,θEL为起始位置换算到EL上的偏移量,ε为扫描该行时的仰角.RA上扫描速度换算到AZ和EL上公式如下:

式中RRA为星表中设定的RA上的扫描速度,RAZ为扫描速度换算到AZ上的速度,REL为扫描速度换算到EL上的速度.

OTF观测程序从观测星表中读取参数,从频率设置文件中读取接收机设置参数.详细观测流程如下:

(1)给FFTS发送源的信息,包括坐标、源名、视向速度、观测频率、参考通道号,FFTS用波束号、源名、观测时间组合作为保存该遍扫描数据文件的文件名创建数据文件.

(2)接收机杜瓦窗口前挡黑体,积分n秒(n从星表中读取),作为校准数据.

(3)打开黑体,天线跟踪参考位置,积分m秒(m从星表中读取),作为参考位置数据.

(4)天线移动到扫描行的起始位置,起始位置为该行扫描起始点加上一个加速区域,赤经(RA)方向扫描为例,起始点计算公式如(1)–(5)式,加速区域在星表中设定,天线跟踪到起始位置后给ACU发送扫描速度指令.扫描速度分解到AZ和EL方向的扫描速度加源的运行速度,扫描速度的分解公式如(6)–(7)式.当天线进入扫描范围,给FFTS发送开始采集数据的指令.指令包含采样间隔ΔT、采样点数、AZ及EL方向扫描速度.FFTS收到指令后给调度计算机返回应答信息,主控机收到应答信息后开始采集天线位置和时间信息,每隔ΔT采集一次,前后两次的值取平均作为该点的值.

(5)当一行扫描结束后把该行所有点的位置和时间信息发送给FFTS,FFTS用(8)式对该行所有ON点测量数据、之前观测的参考点测量数据和黑体测量数据进行校准,然后保存.

(6)天线移动到下一扫描行起始位置,进行下一行扫描,扫描完设定的行(星表中设定的校准间隔换算得到的行数)后重复第2–5步.

观测时OTF扫描采用了隔行扫描的方法,即隔两行扫描一次,一遍分3次扫描,第1次扫描1、4、7···行,第2次扫描2、5、8···行,第3次扫描3、6、9···行,如图2. 这样做的好处在于:(1)在天线扫描速度、采样间隔、观测参考点时间间隔等参数相同的情况下,一次扫描能快速覆盖一个区域.网格化处理后一个位置的谱线是多次结果平均,每次时间相隔较长,一个网格化范围内的点分布在不同时间段,所以能减小望远镜系统性能和天气变化等因素对整个区域分布的影响;(2)隔行扫描的结果在做网格化处理时,一个网格化区域内相邻几行的数据使用的参考点是不同时间观测的,相邻几行的参考点是独立的,对参考点的积分时间可以减少一些.

图2 隔行扫描示意图(左)和观测过程中扫描轨迹图(右)Fig.2 Left:diagram of interlaced scanning,right:scanning track in an observation

3.3 参数优化

OTF观测需要设置很多参数,如扫描区域和加速区域大小、扫描速度、采样率、参考位置的观测间隔、参考位置的积分时间、黑体的观测间隔等.这些参数有一些理论上的范围,但在范围内每个望远镜根据自身特点需要优化,参数优化可以提高观测效率.由于超导成像频谱仪波束间隔为3′左右,如果扫描区域小于6′,9个波束各自的扫描区域没有重叠的部分,所以扫描区域大于6′时9个波束才会有重叠的区域.扫描区域太大时,扫描一遍要很长时间,这不利于望远镜的控制和数据处理,所以扫描一个大区域时需要分为几个小区域,观测后再拼接,扫描区域选10′–30′为宜.为了避免由于采样不够而导致信息的丢失,由Nyquist采样定理可知,扫描行间隔(θRow)不能大于天线的理论半功率波束宽度(HPBW)的一半,13.7 m望远镜在115 GHz的HPBW为50′′,所以θRow要小于25′′,θRow太小则会不利于与扫描方向采样间隔的匹配,并且使得扫描一遍所用的时间太长,通常设为15′′.数据采样率,即单次积分时间(TDUMP)应该大于接收机后端获取一组有效数据所需要的最短时间.需要综合考虑观测时间的长度、信噪比以及扫描区域的大小等因素.为了使扫描方向的采样间隔和θRow一致,扫描速度R用(9)式计算:

TDUMP分别设为0.1 s、0.2 s、0.3 s、0.4 s、0.5 s,选取点源IRC+10216做测试,网格化间隔设为30′′,得到的中心点积分强度和观测时间归一化后噪声如图3,可以看出,采样间隔取0.1 s时强度明显降低,一个原因是采样率为0.1 s时,天线扫描速度为150 arcsec/s,目前位置读取周期是0.05 s,得到的位置误差会较大,导致波束污染.所以不建议使用0.1 s.采样间隔从0.2 s到0.5 s时,强度仍逐渐增大.0.2 s、0.3 s、0.4 s、0.5 s的归一化噪声接近,可以根据扫描区域大小选择合适的采样间隔,使得扫描1遍的时间不要大于2 h.

由于13.7 m望远镜副面采用主动控制[6−7](补偿重力变形)和自适应调焦[8](补偿温度变形),根据仰角和温差进行副面自动调节,如果观测参考位置和观测源时的副面位置不同,会导致谱线基线不平,所以每次在观测黑体时发送命令给副面控制计算机进行调节,之后观测参考位置和对源扫描时副面不动.如果长时间不动副面会影响指向,所以校准间隔受副面控制限制.根据补偿公式计算,在特殊位置(如观测DEC接近37°的源仰角大于80°时),要把指向误差控制在5′′内,两次OFF间隔最大3.5 min,其他情况下允许的时间间隔会比3.5 min长,但是校准时间间隔过长,接收机增益漂移等会影响谱线质量,因此校准时间间隔不宜过长,通常采用4 min.

4 数据预处理相关参数测试及处理过程

数据预处理是在Linux系统的服务器上,用C++语言编写数据修正软件,对每个采样点数据进行位置、频率和效率修正,然后将9个波束的数据合并到一起,生成一个CLASS格式文件,CLASS是GILDAS(http://www.iram.fr/IRAMFR/GILDAS)软件的一部分.然后用CLASS脚本做坏谱线剔除和Gridding处理.在多波束OTF观测中,需要得到每个波束指向天空的位置相对于中心波束在AZ和EL方向的间隔(称作波束间隔).各个波束的效率会不相同,还需要测试各个波束相对于中心波束的效率(称作相对效率),波束间隔和相对效率通常在每年夏季对接收机维护后测试1次.下面介绍测试方法和数据预处理具体过程.

4.1 波束间隔测试

理论设计时每个波束与中心波束的间隔是相等的,由于馈源喇叭是独立安装的,每次安装的实际值跟理论值是不一样的,需要通过测试得到实际的间隔.测试时先要用中心波束做指向校准,然后用木星等信号强的连续谱点源做AZ-EL方向的OTF扫描观测,即测试每个波束的方向图.方向图扫描范围要覆盖所有波束,如超导成像频谱仪的理论间隔是3′,扫描范围设为12′.对每个波束的数据做二维高斯拟合,得到每个波束的中心位置(ΔAZ,ΔEL),计算出各个波束相对于中心波束的间隔(DAZ,DEL).测得的望远镜的方向图如图4,波束间隔如表1.

图3 设置不同TDUMP观测点源IRC+10216得到的中心点积分强度I和观测时间归一化后的噪声rms(root mean square)Fig.3 Integrated intensity I of the central point toward the point source IRC+10216 and the rms(root mean square)of noise normalized by the observation time at different TDUMP

图4 用OTF模式对木星扫描测得的望远镜方向图,颜色棒为Fig.4Continuum beam pattern of telescope measured by means of scanning Jupiter in the OTF mode,the color bar is

表1 测得的超导成像频谱仪波束间隔Table 1 Beam interval of the superconducting imaging spectrometer

4.2 相对效率测试

效率测试只用中心波束做了测试,由于每个波束的效率不一样,需要测试其他波束强度与中心波束强度的比值,在对数据做效率修正时其他波束的效率为中心波束的效率乘以对应的相对效率.测试方法:在得到波束间隔后,用9个波束轮流对标准源(如NGC 2264、S140)做谱线单点位置调制观测,用某个波束观测时在中心波束观测的天线位置AZ和EL上分别加上前面测得的对应波束间隔(DAZ,DEL).用同一轮观测中其他波束观测的强度得到一组相对效率,为了减小误差,需要多个源的多组效率平均,测得的各个波束的效率如表2.

表2 测得的超导成像频谱仪相对效率Table 2 Relative efficiency of the superconducting imaging spectrometer

4.3 数据修正

4.3.1 天线位置修正

由于扫描时天线位置在调度机采集,谱线数据在FFTS频谱仪采集,有时由于两台计算机通讯不畅造成采样不同步,为了解决此问题,每个采样点都标记时间戳.FFTS计算机和调度机每半小时通过Linux命令ntpdate校准计算机时间,保证计算机时间准确.每个采样点记录计算机时间,与调度机记录的时间进行比较,计算两台计算机采样不同步的时间差,如果时间差小于3 s,通过扫描速度修正AZ、EL,如果时间差大于3 s时,该数据要剔除.用波束间隔和每个采样点对应的LST、儒略日、AZ、EL可以计算各个波束每个采样点对应的实际坐标,用实际坐标和中心坐标计算θRA和θDEC,后面在文件中记录中心坐标、θRA和θDEC.

4.3.2 频率修正

在计算本振频率时,第2级变频晶振频率按理论值计算,实际上每个波束的第2级变频晶振存在偏差,通过频率定标的方法确定出每个波束晶振的实际值,计算出晶振频率偏差导致的速度偏差.另外由于OTF观测扫描一个较大的区域,每个位置的频率多普勒修正项中的本地静止参考系速度VLSR会有差别,控制软件难以实时逐点改变本振频率来修正.OTF观测时只能每行设置一次本振频率,计算本振频率用中心点的VLSR0.对各扫描点n根据其LST、RA、DEC计算出对应的VLSRn,然后用VLSRn减VLSR0计算该点的速度修正值.上述两方面的修正可以通过修改CLASS文件头信息中的Velocity值来实现.

4.3.3 效率修正

用相对效率和主波束效率对每个波束的数据进行效率修正,得到与同类望远镜可比的“观测辐射温度”,即.

4.4 单个坏谱线剔除

4.4.1 剔除坏通道

FFTS有16 k个通道,有的通道会有坏数据,如图5.通过对数据的分析,根据13.7 m望远镜观测数据的特征,采用一些判据,把坏通道的数据用相邻的通道数据替换.如某个通道的值大于300 K或小于−300 K,或某个通道比前后的通道都大50 K或都小50 K,认为这个通道是坏通道,用旁边通道的值替换该通道的值.如果一条谱线上坏通道大于10个,该条谱线质量标记为差,CLASS中有数据质量标志,处理时可以设定数据质量标准把质量差的数据剔除.

4.4.2 剔除基线差的谱线

有时由于接收机增益不稳等原因导致谱线基线不平整,如图6,这些谱线需要剔除.由于每条谱线积分时间很短,噪声很大,不好判断基线是否平整,把谱线做32个通道平滑处理(第1−32个通道平均,第2−33个通道平均,依次类推),然后取速度为−500 km·s−1至−200 km·s−1和200 km·s−1至500 km·s−1两段数据做线性拟合,计算rms,当rms大于限定值时标记该点数据质量为差,限定值根据当时系统温度用经验公式计算.

图5 谱线上的坏通道Fig.5 A bad channel in a spectrum

图6 基线差的谱线Fig.6 A spectrum with a bad baseline

4.5 Gridding处理

OTF观测的过程中进行了很密集的空间采样,这些采样点不一定是等间隔分布,需要进一步的处理,即Gridding,其过程是用卷积函数作加权平均,以获取通常形式的格点数据.CLASS软件具有Gridding功能,用CLASS软件实现Gridding过程.首先根据不同频率的实测值设置望远镜波束大小,重采样间隔设为30′′,然后执行table命令,该命令将当前索引中的所有谱线构建成GILDAS表,再用xy map命令对表中数据采用高斯函数做卷积计算,生成光谱Cube数据及相关的权重图像,然后用vector fits命令生成国际通用的FITS Cube文件,各命令用法详见GILDAS网站介绍.

5 实验结果与分析

为了验证观测过程和数据处理过程的正确性,选取结构比较特殊的源S140进行多遍OTF观测.用4遍观测的9个波束的数据分别进行成图,比较9个波束的分子云分布图,如图7,采样间隔为0.3 s,速度积分范围为−13 km·s−1至−2 km·s−1,比较可以看出:各个波束的分子云分布图形状一致,中心点位置正确,中心点谱线强度一致.每个波束的分布图与9个波束合并后的分布图(如图8)一致,最大值出现的位置一致,相同位置的谱线轮廓和强度一致.经过测试,在不同方位多遍OTF观测得到分子云的分布同样具有一致性,且与FCRAO观测结果(如图9)十分吻合[10].用OTF观测Gridding后中心点的谱线与该位置用位置调制方法做单点观测得到的谱线进行比较,二者强度和轮廓一致,峰值强度差值在rms范围内,如图10,表明我们的观测模式和数据处理方法是可靠的.

图7 9个波束获得的S14013CO的分子云分布图Fig.7 13CO integrated intensity maps of S140 by 9 beams

图8 S14013CO积分强度图Fig.8 13CO integrated intensity map of S140

图9 FCRAO S14013CO积分强度图[10]Fig.9FCRAO13CO integrated intensity map of S140[10]

图10 OTF观测中心点谱线与单点观测谱线比较Fig.10 Spectra from the central point in an OTF observation and a single point observation

6 结论和展望

超导成像频谱仪和OTF观测方法的应用,使得13.7 m望远镜观测效率得到很大提升,10′×10′的区域成图用时40 min得到的图像,如果采用原有的单像元系统,覆盖同样的区域,采用逐点扫描方式得到质量接近的结果则需要700 min[4].OTF观测方法从2011年5月用于13.7 m望远镜谱线成图观测,目前主要用于对北天银河系平面±5°范围内的CO同位素3条分子谱线进行大尺度巡天,其他课题成图观测中也得到广泛的应用.目前的数据预处理过程部分采用CLASS软件,先把数据写成CLASS格式,再用CLASS软件读取,多次读写文件,数据预处理花费时间较长,预处理方法需要进一步优化.

致谢OTF观测前期杜福君做了调研工作,系统开发中银河画卷项目组成员给予很多帮助和对观测结果进行测试验证,青海观测站的同事给予大力支持,在此表示真挚的感谢.

[1]Mangum J G,Emerson D T,Greisen E W.A&A,2007,474:679

[2]Sawada T,Ikeda N,Sunada K,et al.PASJ,2008,60:445

[3]董健,吴亚军,袁瑾,等.天文学进展,2016,34:212

[4]杨戟,曹凝.中国科学院院刊,2011,26:478

[5]Shan W L,Yang J,Shi S C,et al.ITTST,2012,2:593

[6]左营喜,李阳,孙继先,等.天文学报,2011,52:152

[7]Zuo Y X,Li Y,Sun J X,et al.ChA&A,2011,35:439

[8]孙继先,左营喜,杨戟,等.天文学报,2014,55:246

[9]Dame T M,Hartmann D,Thaddeus P.ApJ,2001,547:792

[10]Ridge N A,Wilson T L,Megeath S T,et al.AJ,2003,26:286