2. 中国科学院国家天文台, 北京 100101;
3. 中国科学院大学, 北京 100049
2. National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China;
3. University of Chinese Academy of Sciences, Beijing 100049, China
光学波段天文图像通常由望远镜配备的电子耦合器件(Charge Coupled Device, CCD) 芯片相机采集,主要采用FITS (Flexible Image Transport System) 格式进行存储。为了测量图像上特定星像的亮度或视星等,我们一般采用较差测光或图像相减的方式[1-2]。对于前者,先在观测图像(测光图像) 上选取一些参考星,通过比较参考星的仪器星等和这些参考星在巡天星表中的视星等,得到图像零点,再对上述特定源的仪器星等进行零点修正,从而得到该源基于巡天星表的视星等。对于选取的参考星,我们需要通过其在CCD图像中的像素坐标计算出对应的天球坐标,再与巡天星表中的对应目标进行匹配。观测图像的位置定标越准确,图像中星像与参考星表中恒星的坐标匹配结果越可靠,否则可能出现误匹配(图像某星像与匹配上的星表中的恒星不是同一颗星) 的情况,从而影响零点的准确度。同时,获得并报告特定源的天球坐标信息也是测光数据处理的主要目标之一。显然,观测图像上各星像的位置定标是整个测光数据处理过程中非常重要的一环。
从图像像素坐标到天球坐标的转换关系的系数以及相关信息记录为FITS文件头中的一系列特定关键字的数值和备注信息,这些关键字称为图像的世界坐标系[3] (World Coordinate System, WCS) 信息。根据文[4],从图像像素坐标到天球坐标的转换可以分为3步:(1)从像素坐标(Pixel Coordinates) 转换到以度为单位的投影平面坐标(Projection Plane Coordinates);(2)从投影平面坐标转换到本地球面坐标(Native Spherical Coordinates);(3)从本地球面坐标转换到最终的天球坐标(Celestial Spherical Coordinates)。第2步操作视投影类型而定,有固定的转换公式,在图像中常使用心射正切投影(TAN[4])。第3步为不同球面坐标系下的坐标转换,转换公式也是固定的。因此,天测位置定标实际要计算的就是第1步中的转换系数。
理想情况下,从像素坐标到投影平面坐标的转换可以由线性变换近似描述,包括坐标的平移、旋转、伸缩、切变等,但由于大气和仪器等因素(如大气折射[5]、大气湍流、望远镜CCD缺陷、光路缺陷和重力形变等),望远镜拍摄到的图像相比理想的投影平面存在一定程度的扭曲。为了保证图像各个区域坐标具有足够的准确性,一般的方法是引入高阶多项式作为修正项拟合图像扭曲,以修正坐标偏差。若望远镜视场较小且仪器状态良好,则对所拍摄图像定标时仅计算线性转换系数(后文简称线性定标) 就可以满足需求,但对于巡天图像这类视场较大的图像,边缘部分扭曲比较明显,此时就需要额外计算高阶修正系数(后文简称扭曲修正) 以得到理想的定标结果。
当前,Astrometry.net [6]和SCAMP[7]是使用较多的两款天测位置定标软件,前者主要特点为在没有任何关于图像位置信息时可以计算一套完整的世界坐标系信息,而后者主要功能为对一套初始的世界坐标系信息进行修正。本文首先介绍两款软件的定标原理,接着以ZTF (美国的一项正在进行中的时域巡天项目) 巡天[8]的科学图像为测试图像,计算并比对了在输入源表相同的情况下,以4种不同方式运行两款软件得到的天测定标质量。文中使用的Astrometry.net版本为0.89,SCAMP版本为2.10。
1 天测位置定标原理介绍如前所述,从像素坐标到天球坐标转换的关键在于像素平面坐标到投影平面坐标的转换。若不考虑扭曲修正,则两者之间的转换关系由FITS头文件中的CDi_j给定,i和j取值1或2。线性转换公式为
$ \left(\begin{array}{l} x \\ y \end{array}\right)=\left(\begin{array}{ll} C D_{1 \_1} & C D_{1 \_2} \\ C D_{2 \_1} & C D_{2 \_2} \end{array}\right) \times\left(\begin{array}{l} u \\ v \end{array}\right), $ | (1) |
其中,(u, v) 为以参考点为原点的像素坐标;(x, y) 为以同一参考点为原点的投影平面坐标,参考点由关键字CRPIX1和CRPIX2给定。若考虑扭曲修正,依据所使用的扭曲修正类型,先修正像素坐标(u, v) 再进行线性转换,或者先进行线性转换再修正投影平面坐标(x, y)。
天测定标软件一般通过图像星表和参考星表之间的比较分析来计算坐标转换中的各项系数。对于Astrometry.net,使用者既可以输入图像由软件完成找星生成图像星表,也可以直接输入准备好的图像星表由软件直接进行定标。而对于SCAMP,使用者需要提前使用软件SExtractor[9]生成图像星表,再输入SCAMP进行定标修正。
1.1 Astrometry.netAstrometry.net是用C语言编写的开源免费软件,可以在Linux系统和MacOS系统上安装(见https://astrometrynet.readthedocs.io/en/latest/),安装后可以通过solve-field命令运行。该软件通过匹配图像星表中的星标(asterism) 和参考星表中的星标计算图像的世界坐标系信息。一个星标由4颗星组成,根据这4颗星的坐标可以计算出一个代表它们之间相对位置关系的哈希码;如果这4颗星既在图像星表中也在参考星表中,则由它们的像素坐标计算的哈希码和由他们在参考星表中的天球坐标计算的哈希码是接近的(由于存在位置测量误差和扭曲,它们的哈希码不完全相同)。
Astrometry.net根据不同的光学巡天星表构建不同尺度的遍布全天的星标,这些星标的哈希码以及组成相应星标的星的天球坐标存储为一系列index文件(http://data.astrometry.net/) 以供查找比对。Astrometry.net根据亮度排序,依次对图像星表中的星像构建星标并计算哈希码,在index文件中寻找哈希码相近的星标作为匹配星标。Astrometry.net根据找到的匹配星标计算世界坐标系信息,据此计算图像上其他星的天球坐标,并与index文件中其他星的天球坐标进行比较,通过贝叶斯决断方法判断该匹配是否为真。若匹配为真,则完成计算;若为假,则继续验证其他的匹配星标。
作为可选项,Astrometry.net可以根据匹配上的星计算图像扭曲修正系数,采用的扭曲修正多项式为SIP (Simple Imaging Polynomial) [10],修正对象为像素坐标(u, v),公式为
$ \left(\begin{array}{l} x \\ y \end{array}\right)=\left(\begin{array}{ll} C D_{1 \_1} & C D_{1 \_2} \\ C D_{2 \_{1}} & C D_{2 \_{2}} \end{array}\right) \times\left(\begin{array}{l} u+f(u, v) \\ v+g(u, v) \end{array}\right), $ | (2) |
其中,f(u, v)和g(u, v) 分别为两方向坐标的修正多项式函数,
$ f(u, v)=\sum\nolimits_{p, q} A_{\_p\_ q} u^p v^q, p+q \leqslant A_{\_ \text {ORDER }}, $ | (3) |
$ g(u, v)=\sum\nolimits_{p, q} B_{\_p\_ q} u^p v^q, p+q \leqslant B_{\_ \text {ORDER }} \text {, } $ | (4) |
A_ORDER和B_ORDER为两个方向的修正阶数,一般是相等的,可以通过命令行参数设置;A_p_q和B_p_q为多项式中的系数。
1.2 SCAMPSCAMP是用C语言编写的开源免费软件,可以在Linux系统和MacOS系统上安装(见https://github.com/astromatic/scamp/),安装后可以通过scamp命令运行。该软件需要使用者提供一个初始的世界坐标系信息,通过独立的两步完成对世界坐标系信息的修正。(1)模式匹配(Pattern matching),目的是对初始世界坐标系中的线性定标系数进行修正;(2)计算天测解(Astrometric solution),目的是在给定线性定标的基础上计算扭曲修正。使用者可以根据需要,只运行其中一步,或者两步都运行。
在模式匹配中,SCAMP首先依据初始世界坐标系信息,把参考星表中的天球坐标投影转换为像素坐标。通过互相关分析比较图像星表像素坐标对和参考星表像素坐标对的间距和方位角的二维分布(见文[7]第3节),SCAMP可以得到初始世界坐标系中的像素比例尺和旋转角度的修正量,并据此更新世界坐标系和参考星表像素坐标。通过同样的方法比较图像星表像素坐标和参考星表像素坐标在像素平面上的分布,SCAMP可以得到两套坐标之间的平移偏差,从而对初始世界坐标系中的指向信息进行修正。
在计算天测解中,SCAMP对图像星表和参考星表的源进行匹配,通过最小化加权坐标差值平方和(见文[7]第4节) 的方法计算扭曲修正系数。SCAMP采用的扭曲修正多项式为TPV[5],修正阶数最高可设为7阶,修正对象为投影平面坐标(x, y),修正公式见文[5]中的(1)式和(2) 式。
2 定标结果对比实验 2.1 测试图像选择本文随机选择了100个均匀分布于北天的坐标,在每个坐标位置处选择一张具有较好视宁度和较深极限星等的r波段ZTF巡天科学图像作为测试图像(下载方式见https://irsa.ipac.caltech.edu/docs/program_interface/ztf_api.html)。每张图像的视场约为52′ × 52′,像素尺寸为3 072 × 3 080,像素比例尺为1.012角秒/像素。图 1展示了测试图像的半高全宽、极限星等(与参考星表较差后对应波段的视星等) 以及用于天测定标的源(见2.2节) 数量分布直方图。半高全宽和极限星等均从图像的FITS头文件中获取。
每张ZTF科学图像的FITS头文件中带有准确的世界坐标系信息。为了量化测试图像的扭曲程度,我们利用图像的世界坐标系信息,计算每个像素坐标考虑扭曲修正和不考虑扭曲修正分别对应的天球坐标之间的偏差,取所有像素位置的偏差距离的最大值作为图像扭曲程度的量化指标。图 2展示了随机选取的一张图像的扭曲情况。在此指标下,所有测试图像中扭曲程度最大为1.46″ (约1.44像素),平均为0.75″ (约0.74像素)。图 3展示了图像扭曲程度的分布直方图。
2.2 定标流程和参数设置首先删除图像已有的世界坐标系信息,使用SExtractor对图像进行找源,并基于该软件给出的每个源的测量结果筛选适合用于定标的源。找源时的重要参数设置见表 1。筛选时,对源的半高全宽进行3倍标准差筛选以剔除半高全宽过大的源,仅保留FLAGS为0的源以去掉测量结果可能不准确的源,同时去掉椭率大于0.5、信噪比低于10的源,最后得到用于天测定标的图像星表。
Parameter | Value | Explanation |
DETECT_MINAREA | 5 | The minimum area above the detection threshold is required |
DETECT_THRESH | 5.0 | Detection target threshold |
ANALYSIS_THRESH | 5.0 | Analyze the target threshold |
FILTER | N | Use filter or not |
PHOT_AUTOPARAMS | 2.5, 3.5 | MAG AUTO parameter, Kron factor and minimum radius |
依次通过4种不同流程进行位置定标,分别为:
(1) 运行solve-field命令,一次性进行线性定标和扭曲修正;
(2) 首先运行solve-field进行线性定标(使用-T参数),再次运行solve-field,在第1次运行的基础上(使用-v参数) 进行包括扭曲修正在内的计算;
(3) 首先运行solve-field进行线性定标,接着运行SCAMP直接进行扭曲修正;
(4) 首先运行solve-field进行线性定标,接着运行SCAMP进行线性定标的修正,最后运行SCAMP进行扭曲修正。
计算的扭曲修正均为4阶,与图像自带的修正阶数一致。这4种流程按次序分别编号为流程1到流程4。不同流程中运行solve-field时固定的命令行参数选取和设置见表 2。不固定的参数只有搜索中心坐标的设置、是否计算扭曲修正及其阶数的设置和是否验证已有世界坐标系的设置。不同流程中运行SCAMP时固定的重要参数设置见表 3,不固定的只有是否进行模式匹配和扭曲修正。Astrometry.net的参考index文件选为基于Tycho-2星表和Gaia DR2星表的5200系列,SCAMP使用的参考星表为Gaia DR2星表。所有计算过程在内存为16 GB、芯片为Apple M1 Pro的14英寸MacBook Pro笔记本电脑上运行。
Command line parameter | Value | Explanation |
--scale-low | 0.7 | Lower limit of image field of view |
--scale-high | 0.9 | Upper limit of image field of view |
--scale-units | degwidth | The unit of the upper and lower limits of the field of view |
--radius | 0.8 | Search radius |
--crpix-center | --1 | Specify the center of the image as the reference point |
--width | 3 072 | Image width (pixel) |
--height | 3 080 | Image height (pixel) |
--x-column | XWIN_IMAGE | Enter the column name for the horizontal coordinate in the image star list |
--y-column | YWIN_IMAGE | Enter the column name for the vertical coordinate in the image star list |
--sort-column | FLUX_AUTO | The name of the column to sort the image star list |
1:该命令行参数不用取值,选定即可。 |
Parameter | Value | Explanation |
ASTREFMAG_LIMITS | 12.0, 18.0 | Refer to the magnitude range used for scaling in the star list |
PIXSCALE_MAXERR | 1.2 | Maximum error multiple of pixel scale |
POSANGLE_MAXERR | 5.0 | Maximum azimuth deviation (degree) |
POSITION_MAXERR | 1.0 | Maximum position deviation (Angle minute) |
PROJECTION_TYPE | SAME | Projection type |
DISTORT_KEYS | XWIN_IMAGE, YWIN_IMAGE | The column name in the image star list on which the distortion correction depends |
DISTORT_DEGREES | 4 | Distortion correction order |
NTHREADS | 1 | Number of running threads |
考虑到不同使用环境对天测定标的准确性要求不同,本文使用定性和定量两个指标衡量天测定标的质量。定性指标选为定标后图像星表源坐标是否在全图范围内与参考星表源坐标相匹配,匹配阈值选为0.5倍图像半高全宽,用于匹配的参考星表选为Gaia DR2星表。定量指标则选为定标后的图像源坐标与参考星表中对应源坐标的均方根偏差。具体计算方法如下。
(1) 使用SExtractor对定标后图像分别以高低阈值进行两次找源,高探测阈值和分析阈值选5.0,低探测阈值和分析阈值选1.0。以与2.2节中相同的筛选条件对高阈值星表进行筛选,再与参考星表进行坐标匹配,接着将参考星表与低阈值图像星表匹配。如果与高阈值星表源配对的参考星表源在低阈值星表中配对的源不是当前高阈值星表源,则判定为误匹配。在高阈值星表中去掉所有误匹配的源,计算余下源坐标与匹配上的参考星表源坐标之间的距离,作为图像星表源坐标的偏差。需要指出的是,SExtractor不能识别由SIP多项式描述的扭曲修正关键字,因此,对于Astrometry.net给出的定标结果,我们先使用python包sip_tpv(https://github.com/stargaser/sip_tpv) 将SIP关键字转换为TPV关键字,再使用SExtractor进行找源。
(2) 对图像划分为10行10列的网格阵列,将每个网格内源的坐标偏差除以0.5倍图像半高全宽,再剔除5倍标准差(标准差选为84.13百分位与15.86百分位差值的一半) 之外的异常值后,取平均值作为当前网格内源坐标准确性的衡量,最后得到的数组即可反映图像各区域坐标的偏差情况。图 4展示了一张图像只进行线性定标和扭曲修正后的坐标偏差图。如果数组内所有值小于1,则认为图像星表在全图范围内与参考星表相匹配,定性为定标质量好;如果数组内大于1的值的个数大于0且不超过10,则定性为定标质量中等;否则定性为定标质量差。
(3) 计算剔除异常值后所有源的赤经赤纬坐标偏差(定标后的坐标减去参考星表的坐标) 以及相应的均方根,作为该图像天测定标质量的定量指标。图 5展示了所有图像上的源只进行线性定标和进行扭曲修正后的坐标偏差分布图。
除了以上提到的指标,赤经赤纬平均偏差及偏差的标准差与星等的关系图也能直观地展示定标质量的好坏。图 6展示了所有图像上的源在扭曲修正后赤经赤纬偏差随星等的分布图,包括各星等处的偏差中值及1倍标准差范围。由图 6可见,两方向的坐标偏差在各星等处较均匀地分布在零点左右。
3 定标结果比对本节从总计算耗时、定性的星表匹配情况和定量的坐标均方根偏差3个方面对不同定标方式得到的结果进行比较。比较对象除了2.2节中提及的4种定标流程的结果,也包括Astrometry.net仅线性定标的结果和SCAMP仅线性定标修正的结果,以及ZTF图像已有的定标结果。表 4记录了这些定标结果的平均计算耗时,定标质量分别为好、中等、差的图像占比,以及100张图像平均赤经赤纬方向的均方根偏差。
Calibration process | Average time cost/s | Fraction of good quality/(%) | Fraction of medium quality/(%) | Fraction of bad quality/(%) | Average RMS along R.A./mas | Average RMS along Dec/mas |
linear calibration by Astrometry.net | 0.72 | 29 | 40 | 31 | 329.18 | 314.60 |
improved linear calibration by SCAMP | 4.59 | 68 | 32 | 0 | 184.17 | 169.82 |
Process 1 | 1.64 | 83 | 8 | 9 | 177.81 | 173.26 |
Process 2 | 5.95 | 96 | 4 | 0 | 73.39 | 69.99 |
Process 3 | 1.00 | 100 | 0 | 0 | 68.68 | 64.36 |
Process 4 | 5.59 | 100 | 0 | 0 | 67.80 | 63.69 |
ZTF calibration | — | 100 | 0 | 0 | 68.66 | 64.58 |
定标流程4中包含SCAMP对Astrometry.net给出的线性定标结果进行修正的步骤,因此可以对修正前后的线性定标质量进行比较。图 7展示了修正前后定标质量分别为好、中等和差的图像数量分布。图 7表明,经过SCAMP修正后,定标质量好的图像占比由29%升高至68%,而质量差的图像占比由31%降至0%。图 8展示了修正前后赤经赤纬方向均方根偏差的分布。由图 8可以看出,经过修正后,两方向均方根偏差分布更集中且均值更小。图 9展示了运行solve-field进行线性定标和单线程运行SCAMP进行线性定标修正的耗时分布。由图 9可以看出,Astrometry.net一般在1 s内即可完成图像的线性定标,而SCAMP则需要4~5 s才能完成线性定标的修正。
3.2 4种流程定标结果对比我们对4种不同定标流程所得结果的质量进行对比。图 10展示了4种流程定标质量分别为好、中等和差的图像数量分布。由图 10可以看出,流程1的定标结果中质量好的占比最低,仅为83%;流程2的定标结果中质量好的占比提高至96%,但仍没有达到100%。相比之下,由SCAMP参与定标的流程3和流程4所得定标结果中质量好的占比均达到了100%。这表明搭配使用Astrometry.net和SCAMP进行定标比仅使用Astrometry.net定标的效果要好。图 11展示了4种流程定标后赤经赤纬方向均方根偏差的分布。由图 11可以看出,流程2,3和4所得定标结果的均方根偏差分布相差不大,集中程度和平均偏差均优于流程1的结果。值得注意的是,流程3和4定标结果的相差很小,说明在使用SCAMP做扭曲修正前,是否对Astrometry.net给出的线性定标结果进行修正对最终定标结果的准确度影响较小。图 12展示了4种定标流程的总耗时分布。图 12表明,流程3的总耗时最少,大部分图像在1 s可以完成定标,这说明SCAMP计算扭曲修正耗时极短;流程2的总耗时最长,一般需要5~6 s内才能完成图像的定标,这是因为在流程2中,Astrometry.net在第2次运行的过程中使用多个不同尺度的index文件进行验证并计算扭曲修正,因此耗时较长。
4 讨论结合3.1节和3.2节的对比结果以及两款软件定标原理可知,Astrometry.net相比于SCAMP的主要优点为在没有预先进行世界坐标系估计的情况下以较快速度对图像进行线性定标,缺点在于若要进行扭曲修正,一次性定标虽然耗时短,但质量较差,分两次定标则虽然质量较好,但耗时较长;与Astrometry.net相比,SCAMP的主要优点在于能够快速基于给定的线性定标结果进行扭曲修正,而对给定线性定标进行修正虽然能够显著提升定标质量,但耗时较长。搭配使用Astrometry.net和SCAMP可以实现快速且准确的天测位置定标。
同时需要指出的是,(1) 如果实际拍摄图像的扭曲程度比本文中测试使用的ZTF图像的扭曲程度严重,那么有可能导致Astrometry.net线性定标质量太差,从而使得SCAMP无法通过直接扭曲修正获得足够准确的定标,此时很可能需要先进行一次线性定标的修正再进行扭曲修正;(2) 本次实验所选择的测试图像均为质量较好的图像,对于深度较浅、视宁度较差的图像(图像的科学价值会随之降低),不同定标方式之间的对比结果可能发生变化;(3)本文在使用Astrometry.net的solve-field命令时输入的搜索中心坐标为根据图像已有世界坐标系信息计算的准确中心坐标,而实际数据处理中往往只能知道图像大致的指向坐标,因此实际的解算耗时可能更长。
5 总结与展望本文介绍了Astrometry.net和SCAMP两款软件的天测位置定标原理,以ZTF的100张巡天科学图像为测试图像,从计算耗时、图像星表与参考星表匹配情况、坐标均方根偏差3个方面比较了使用两款软件以4种不同的定标流程进行定标的结果。研究发现,虽然Astrometry.net能够在没有预先世界坐标系信息估计的情况下快速地一次性进行线性定标和扭曲修正,但定标准确性逊于耗时较长的分两步计算;而采用Astrometry.net计算线性定标搭配SCAMP进行扭曲修正的方法,可以做到在耗时和准确性上均优于仅使用Astrometry.net进行定标。在后续的研究工作中,我们可以进一步探究是否存在更好的选源方式、更好的软件参数设置以进一步提高定标的速度和准确性。
致谢: 感谢中国科学院国家天文台郑捷博士提供关于量化图像扭曲程度的建议,感谢华中科技大学博士研究生朱子佩在论文修改方面提供的建议。
[1] |
郑捷, 赵刚, 王炜, 等. SAGE测光巡天数据处理方法研究[J]. 天文研究与技术, 2019, 16(1): 93–106 ZHENG J, ZHAO G, WANG W, et al. Research on the data reduction of the SAGE photometric survey[J]. Astronomical Research & Technology, 2019, 16(1): 93–106. |
[2] |
潘祖恒, 彭青玉, 陆星, 等. 利用图像相关性搜寻运动目标[J]. 天文研究与技术, 2023, 20(1): 41–49 PAN Z H, PENG Q Y, LU X, et al. Searching moving objects by image correlation[J]. Astronomical Research & Technology, 2023, 20(1): 41–49. |
[3] | GREISEN E W, CALABRETTA M R. Representations of world coordinates in FITS[J]. Astronomy & Astrophysics, 2002, 395(3): 1061–1075. |
[4] | CALABRETTA M R, GREISEN E W. Representations of celestial coordinates in FITS[J]. Astronomy & Astrophysics, 2002, 395(3): 1077–1122. |
[5] | SHUPE D L, LAHER R R, STORRIE-LOMBARDI L, et al. More flexibility in representing geometric distortion in astronomical images[C] // Proceedings of SPIE. 2012. |
[6] | LANG D, HOGG D W, MIERLE K, et al. Astrometry.net : blind astrometric calibration of arbitrary astronomical images[J]. The Astronomical Journal, 2010, 139(5): 1782–1800. DOI: 10.1088/0004-6256/139/5/1782 |
[7] | BERTIN E. Automatic astrometric and photometric calibration with SCAMP[C] // Proceedings of the Astronomical Data Analysis Software and Systems XV ASP Conference Series. 2006. |
[8] | MASCI F J, LAHER R R, RUSHOLME B, et al. The Zwicky Transient Facility: data processing, products, and archive[J]. Publications of the Astronomical Society of the Pacific, 2019, 131: 018003. DOI: 10.1088/1538-3873/aae8ac |
[9] | BERTIN E, ARNOUTS S. SExtractor: software for source extraction[J]. Astronomy & Astrophysics Supplement, 1996, 117: 393–404. |
[10] | SHUPE D L, MOSHIR M, LI J, et al. The SIP convention for representing distortion in FITS image headers[C] // Proceedings of the Astronomical Data Analysis Software and Systems XIV ASP Conference Series. 2005. |