光变曲线自动抽取程序

郑 捷，江林巧

(1. 中国科学院国家天文台，北京 100101；2. 乐山师范学院数理学院，四川 乐山 614000)

随着天文观测技术的发展，时域天文观测研究已经成为天体物理学重要的前沿领域，即通过多波段、时域监测来研究宇宙中各类天体的特征和变化，寻找和探索各类新天体、新现象。在时域天文观测的各类目标中，变星是其中一类重要的观测目标，变星种类繁多[1]，观测研究历史悠久，对天文学的发展有着强大的推动作用。

兴隆观测站隶属于中国科学院国家天文台和中国科学院光学天文重点实验室，是亚洲规模最大的光学天文观测基地。兴隆观测站拥有中小口径光学天文望远镜十余台，比如85 cm望远镜、80 cm望远镜、60 cm望远镜等等。这些望远镜每年面向全球天文工作者开放申请，承担了大量变星的观测研究工作[2]。

在对变星的测光观测研究中，对原始数据进行测光处理，获取多波段、高精度的变星光变曲线是一项基础且重要的工作。多年来，天文学家开发了相应的测光数据处理软件包，比如虽然已经停止维护但仍然广泛使用的IRAF[3]软件中自带各种非常成熟的测光软件包。传统的测光数据处理方式需要手动对原始数据进行逐步操作，数据处理效率低。近年来虽然各类光变曲线处理程序不断涌现[4-6]，自动化程度、易用性等比以前有较大幅度提高，但由于不同的望远镜对观测数据文件的命名以及头文件中各类关键词信息等都未统一，导致很多光变曲线处理程序的通用性欠佳，即程序只适用于某一特定望远镜产生的数据，当应用于其他望远镜时，程序常常需要手动更改参数方可继续运行，甚至无法运行。兴隆观测站作为亚洲最大的光学天文观测基地，多台米级口径望远镜每天都在进行大量变星测光观测，但至今没有完全针对兴隆站米级口径望远镜的测光数据自动化处理软件。

鉴于此，本文针对兴隆米级口径望远镜研发了一套光变曲线自动抽取程序，只要提供原始观测数据的保存目录，程序即可实现数据类型识别、数据预处理、测光、图像对齐、参考星选取、较差定标等步骤，并自动进行光变曲线绘制。该程序命名为Quick Light-Curve Pipeline，简称QLCP。该程序已经发布在国内代码托管平台 “码云(Gitee)”。通过运行程序，观测者可以在对目标进行测光观测的同时，快速对原始数据进行测光处理，获得目标的光变曲线，以便能够准确合理地制定下一步的观测计划，具有一定的科研辅助价值和应用价值。

1 程序整体框架

1.1 支撑软件和环境

QLCP程序采用Python3语言编写，在找源和测光部分调用文[7]开发的天文软件Source-Extractor(以下简称SEx)完成。此外，程序需要Python天文软件包astropy[8]和pystronomy，以及Python通用软件包numpy, scipy, matplotlib, os, sys和re等的支持。由于SEx仅在Linux或macOS操作系统上工作，程序也受同样限制。

1.2 程序结构和总体流程

程序总体结构如图1，分为单图像处理和多图像处理两部分。单图像处理又分为文件类型识别和列表生成、本底合并、平场合并、图像改正、找源和测光5个模块。多图像处理分为图像对齐匹配、参考星选择、较差定标、光变曲线绘制4个模块。此外，程序还包括1个统一调用模块和若干支撑模块。各功能模块之间主要通过输入、输出文件进行耦合。考虑到不同科研课题对数据处理的要求不同，程序的每个模块均可以使用自定义模块替代。

1.3 优 势

程序具备以下优势：

图1 程序模块示意图。(a)单图像处理；(b)多图像处理

(1)自动识别文件。对兴隆站望远镜的观测数据，程序只需要给定原始数据目录，即可自动识别图像类型，分类进行处理。对于兴隆站之外的其他望远镜数据，也只需要简单配置。

(2)图像自动对齐。通过星像匹配算法，高效率、精准实现连续观测图像的对齐工作。

(3)自动选择参考星。如果用户未指定参考星，程序能够自动从识别的源中选择参考星，通常为图像中信噪比较高且在观测期间光度稳定的源。

(4)使用方便。只需要指定原始数据保存路径、处理结果保存路径和目标在图像中的位置，程序即可自动完成光变曲线抽取，也可以进行自定义处理。

(5)模块化与替代性好。程序各步骤之间通过输出、输入文件进行耦合，任一模块均可以通过自定义操作进行替换，可以局部使用该程序。

(6)自定义性强。程序可以通过配置文件进行参数自定义，无需修改源代码或编译即可使用自定义参数对数据进行处理，可以根据不同的科研需求进行调整，方便应用于不同场合。

2 关键问题和解决方法

程序的处理步骤均参考常用的天文观测数据处理方式[9]，并考虑程序适用的设备类型做针对性的设计。找源和测光部分调用SEx完成，其他部分均为自行编程实现。其中astropy, numpy和scipy等软件包提供基础功能支持，例如读写FITS(Flexible Image Transport System)文件、数组计算、统计等。以下对部分关键步骤加以说明。

2.1 文件识别

程序识别观测数据的信息主要依据文件名提供的信息，例如观测类型(本底、平场、科学目标)、观测目标、观测波段和观测序号等，并自动生成文件列表，供后续使用。由于不同望远镜文件名设置不同，因此，程序通过在配置文件中设置正则表达式进行解析。默认配置文件已经针对兴隆各望远镜进行设置。

如果观测数据文件名中未体现观测信息的数据，可以通过人工生成文件列表的方式进行处理。人工生成列表还可以实现对数据的筛选，剔除观测错误的图像。

2.2 预处理

识别并生成文件列表后，程序进行本底和平场合并，并对数据进行改正。本底和平场合并均采用中值，而平场在合并之前根据每一幅平场的中值进行归一化。数据在处理之后，会在指定的输出目录中分类保存。

由于在许多望远镜探测器中并未包括过扫描或者预扫描区域，因此程序并未对此进行处理。此外，目前大部分专业天文望远镜使用CCD或者COMS作为探测器，并通过制冷装置对探测器进行降温，暗流往往处于很低的水平，不需要特意加以改正，因此程序也未进行暗流校正。最后，由于程序针对多幅测光图像进行处理，宇宙线、运动物体轨迹等干扰因素带来的假源在多幅图像识别中会自动剔除，程序也未对此进行改正。

在预处理时，程序会进行FITS文件头重建，根据目标信息、观测时间等，计算JD，HJD以及望远镜观测俯仰角、方位角、月相和月球距离等信息。

2.3 测 光

在图像中找源和测光的工作，由程序调用成熟的测光软件SEx完成，需要用户自行安装。SEx需要的配置文件，在程序中内置了default.sex文件，可以满足大部分情况的需要。如对找源、测光等步骤有特殊配置要求的，用户可以自定义。

SEx输出字段在程序内置的default.param文件中指定，也可以自定义，但至少要求包括目标的中心坐标、仪器星等和误差等信息，星像的半高全宽(FWHM)和伸长率(Elongation)等信息为可选输出。仪器星等默认选用MAG_AUTO字段，也可以自定义。

2.4 图像对齐匹配

处理连续多幅测光图像时，程序先自动实现图像之间的匹配，从而获取每颗星在不同图像上的位置。程序在基准图和待匹配图中分别选出最亮的若干颗星，计算两者之间所有星对的x,y坐标距离。由于星像是随机分布的，坐标距离的众值必然是两图之间的实际位移。得到粗略的x,y位移量之后，以此为基础进一步目标匹配，得到更加精确的位移量。

基准图像默认为文件列表的第1幅图像，也可以指定列表中任意一幅，甚至可以指定非当晚观测数据作为基准。

2.5 自动选择参考星

对于大部分被观测的变星，观测者往往在其周围指定若干颗参考星和检验星。程序允许用户手工指定参考星，但是作为一个自动数据处理软件，程序还可以自动从连续观测的多幅图像中找到亮度稳定的亮星作为参考星。具体做法是根据图像匹配结果，计算每一幅图像中仪器星等和基准图像的仪器星等的差值，选择其中信噪比较高的星，并剔除异常点之后计算均值，作为该幅图像和基准图像之间的零点差。随后将图像仪器星等改正到基准图像的同一水平。最后计算这些高信噪比星改正后星等在各图之间的标准差，选择标准差最小的星作为参考星。

自动选择参考星能够大幅度提高数据处理效率，但是也存在缺陷。参考星的选择只根据当晚的观测数据得出，有可能误选长周期变星，但是这对于生成当天的光变曲线没有影响。

3 程序使用和性能

3.1 程序发布

程序已经发布在国内代码托管平台 “码云(Gitee)”上，地址为https://gitee.com/drjiezheng/qlcp/，其中包括程序代码、使用说明等。

3.2 使用方式示例

程序下载后可以放在当前目录下使用。本文以兴隆站60 cm望远镜在2022年5月14日对食双星UY UMa的V波段观测数据进行处理为例，说明处理过程，原始数据保存于目录/data/raw/20220514_60/，处理结果保存于/data/raw/20220514_60/目录中，并且通过人工检查第1幅图像，确认目标星中心点大约在(995, 1 000)，当晚对该目标进行了8.2 h(0.342天)观测，共获得484幅V波段数据，此外还有10幅本底和5幅晨昏光平场。

在配置文件中，除数据处理选项外，程序对观测台站和观测设备也进行了定义。根据兴隆站60 cm望远镜的文件命名规则，观测文件命名方式为UYUMa_V_001.fit，用下划线分割目标名称、波段、观测序号等信息，因此相应的配置文件包括：

site_lon=117.57722 # 117.34.38 #观测站的经度

site_lat=40.395833 # +40.23.45 #观测站的纬度

site_ele=960 #观测站的高度

filename_temp="(?P[^-_]*)_(?P[a-zA-Z]{0,1})_{0,1}(?P[0-9]{3}).fit"。

如果将程序应用于其他台站，只需要修改相应台站信息即可。根据数据文件名命名规则修改filename_temp关键字指定的正则表达式，例如对于兴隆站85 cm望远镜，文件名命名方式为20220514UYUMa-0001V.fit，对应正则表达式为"(?P[0-9]{8})(?P[^-_]*)-(?P[0-9]{4})(?P[a-zA-Z]{0,1}).fit"。

以上述60 cm望远镜观测数据为例，最简单的调用方式如下：

import qlcp22# 引入程序包qlcp22.do_all(# 该程序的统一调用入口 ini_file="xl60cm.ini", # 指定目标望远镜配置文件 raw_dir="/data/raw/20220514_60/",# 指定原始数据目录 red_dir="/data/red/20220514_60/",# 指定输出文件目录 starxy=[(995, 1000)],# 指定目标星的位置)

先引入包，随后调用do_all函数自动完成全部步骤。函数参数最少为以上4个，分别是配置文件(ini_file)、原始数据文件路径(raw_dir)、输出结果路径(red_dir)和目标星位置(starxy)，其余参数采用默认配置。 程序完成从识别文件到最终生成光变曲线的全过程。程序可选的参数还包括对处理步骤的选择、参考星的指定等，以及对处理步骤的细节进行限制的参数。整个处理过程约5 min，最终结果见下一节介绍。

3.3 输出结果

程序各个模块之间通过输出文件和输入文件进行耦合，各模块输出文件主要是处理的中间结果。主要输出文件见表1。

3.4 最终输出数据

图2和图3是上述数据进行自动处理的结果，其中图2标注了目标星(T-0)和自动选择的参考星(T-1到T-5)，图3为得到的目标光变曲线，以及参考星的光变曲线。由于当天观测时气象条件较差，因此得到的光变曲线信噪比不高。从图中可见，当天的参考星光变标准差在0.005到0.010星等之间。

表1 程序的输出文件Table 1 Output files of this program

图2 2022年5月14日利用兴隆观测站60 cm望远镜观测食双星UY UMa的V波段图像，图中标注的T-0为目标星，T-1至T-5为自动选择的参考星

图3 从2022年5月14日利用兴隆观测站60 cm望远镜观测食双星UY UMa的V波段数据中抽取出的光变曲线，其中红色方块为目标星光变曲线，星号为参考星光变曲线，编号顺序与图2中所示一致

3.5 存在的问题

程序只是后续综合程序的基本雏形，针对兴隆观测站米级望远镜的变星观测数据处理进行编制，通过简单配置也可以适用于其他望远镜的变星观测数据处理。在考虑特定使用场景的情况下，程序的功能有一定的局限性，将在未来版本中加以改进。

首先，程序只考虑各幅观测图像之间的平移而不考虑视场旋转，因此该程序不适用于发生场旋的观测。此外，程序并未进行天测定标，因为对于目标明确的变星观测，对星图中目标的具体赤经和赤纬并没有特定的要求，天测等问题将另行解决。由于变星研究通常对目标的视星等没有特别要求，因此该程序也未进行流量定标。实践中较差定标精度实际上往往高于绝对流量定标。最后，由于程序使用了广泛应用的测光软件SEx进行找源和测光，能力和精度受SEx的限制。

4 总 结

本文介绍了从变星连续测光观测数据中抽取光变曲线的程序QLCP，主要包括程序的结构、简单使用方法和处理结果等。该程序可以大幅度简化实测数据的处理，并且具备模块化、可定制等特点。但是该程序作为一个针对特定使用场景程序，有一定的局限性和不足，将在后续的研发中继续改善。

致谢：感谢兴隆观测站工作人员的支持。感谢国家天文台王汇娟、王佳琪、罗常青，云南天文台戴智斌等人在试用该程序时提出的意见和建议。