适用于大视场的快速星像匹配算法

潘祖恒，方雪清，郭碧峰，陆 星，3，彭青玉*

(1. 暨南大学计算机科学系，广东 广州 510632；2. 暨南大学中法天体测量、动力学与空间科学联合实验室，广东 广州 510632；3. 暨南大学物理学系，广东 广州 510632)

在天文观测中，为了方便资料的归算和处理，通常需要将图像中星像的位置和光度等信息与参考列表中对应的信息进行匹配，简称星像匹配。常见的星像匹配方法有基于向量的方法[1]，基于径向和环向特征的方法[2]，基于三角形的方法[3]和基于局部敏感哈希的方法[4]等。比较成功的应用是Astrometry.net[5]，它使用四边形代替三角形，并通过贝叶斯决策准则来选择正确的图像变换。文[6]研究发现，基于向量的方法比基于径向和环向特征的方法更优，具体表现为基于向量的方法耗时更少，匹配率更高。虽然一些研究指出，四边形匹配算法比三角形匹配算法具有更好的效果[4]。但是具有简单结构的三角形仍是较为常用的一种星像匹配算法[1，3，7]。文[8]注意到，基于三角形的方法是以上3种方法[1-3]中匹配率最高，即最稳定的一种方法。但是，该方法列表中的星所构成的三角形数目巨大，造成存储和检索的困难。不仅如此，匹配需要构造全等三角形，匹配全等三角形需要预先知道图像的比例尺。图像的比例尺不是一成不变的，受到各种因素的影响，如焦距、温度、图像扭曲和镜筒弯曲等。

受Astroalign[9]图像对齐算法的启发，本文提出一种基于k-d树和k-means聚类算法的星像匹配算法。该算法利用k-d树进行三角形盲匹配以找到整幅图像区域的整体变换，然后通过k-means聚类算法将整幅图像分割为多个小区域，从而求解小区域范围内的局部变换。如果不存在这样的局部变换，则用整体变换代替。其中，使用k-d树的最近邻搜索算法可以提高匹配的速度，使用k-means聚类算法可以提高匹配的精度。该算法可以间接计算图像的比例尺，从而自适应地体现图像比例尺的微小变化。

1 方 法

因为前人的方法是基于全等三角形的匹配[1]，所以需要预先知道图像的比例尺。虽然获得图像的比例尺并不困难，但是图像的比例尺会因各种因素的影响而发生微小的变化。因此，理论上，星像匹配应使用相似三角形而非全等三角形。为了减小图像比例尺变化和生成大量三角形的影响，本文3次使用k-d树的最近邻搜索算法：第1次使用k-d树生成三角形不变性元组；第2次使用k-d树匹配三角形不变性元组，从而匹配相似三角形；第3次使用k-d树进行最终验证。本质上，k-d树是一种平衡二叉树，也是一种空间划分树。它使用二分的方式将整个k维向量空间不断划分为若干个局部空间。该算法能够在每次搜索过程中进行分支判断，选择满足要求的局部子空间。该方法有效避免了全局空间搜索，以达到快速匹配的目的。为提高匹配的精度，我们使用k-means聚类算法进行区域划分，以达到局部精准匹配的目的，使用心射投影(中心投影)进行天球坐标系到标准坐标系的投影。在不考虑投影误差的前提下，本文假定参考列表的天球坐标投影到了标准坐标，同时也知道了图像中星像的像素坐标(参考星的数目与图像中星像的数目不一定相等)。我们的主要任务是找出共同对象中参考星的标准坐标和图像中星像像素坐标的一一对应关系。

1.1 创建星像列表和读取参考列表

创建星像列表需要对天文图像进行搜星和测光。本文使用DAOFIND[10]程序进行搜星与测光，然后根据星像的亮度进行降序排序，获得星像列表I。列表I中至少包含像素坐标和光通量两种信息。根据天文图像头文件中的视场中心、图像大小和焦距等信息推断图像对应的天球坐标的视场。从Gaia DR3[11]中下载对应的视场中的星像，并按照星等从小到大排序获得参考列表R，列表R的视场远远大于列表I。参考列表R至少包含天球坐标(赤经和赤纬)、自行(实际计算中可以依据观测历元和自行计算出星表在观测历元的近似位置)和星等3种信息。

1.2 生成不变性元组

本文使用相似三角形进行三角形盲匹配。一对匹配好的相似三角形可以求解一个变换。由于匹配的是相似三角形，所以本文算法还可以间接地求解CCD图像的比例尺。为对相似三角形进行盲匹配，本文使用Astroalign算法中的三角形不变性元组。假设三角形的三条边分别为L0，L1和L2，三角形不变性元组可以表示为

(1)

(2)

为了分析三角形不变性元组的特性，直观地了解三角形不变性元组，我们使用天文图像和参考列表中的星像构建三角形来生成不变性元组，并绘制如图1的三角形不变性元组特征图。生成三角形不变性元组的流程如图2。其中，两次使用了k-d树，第1次使用k-d树的目的是找出5个在位置上最近的点，第2次使用k-d树的目的是通过不变性元组找出可能匹配的三角形。选用5个位置最近的点为一组的分组，可以减少三角形生成的数量，且有利于局部匹配。因为横跨小片区域的三角形和横跨大片区域的三角形即使相似，也不可能导出一个正确的变换，而且横跨大片区域的三角形还可能存在图像扭曲。第2次使用k-d树可以快速找出匹配的相似三角形，从而尝试求解一个变换。

图1三角形不变性元组特征图(该图参考文[9])Fig.1 Triangle invariant tuple feature (The figure is referenced from reference [9])

图2 生成三角形不变性元组流程图Fig.2 Flowchart for generating triangle invariant tuple

1.3 导出变换

由图1可知，五角星(☆)和圆圈(○)距离越近，表示两个三角形相似度越高。两个相似三角形的顶点对应关系可以尝试求解一组变换关系。本文使用Astroalign中的RANSAC[12]算法，主要体现为如果星像列表I中的星像进行该变换后，部分三角形能够满足全等匹配关系，说明该变换是一个星像列表I转换到参考列表R的正确变换；否则，重新选择两个相似三角形求解另一组变换关系。重复以上步骤，直到获得一组正确的变换关系。当找到一组正确的变换关系后，即可找出亮星在像素坐标和标准坐标之间的对应关系，最后使用最小二乘拟合出一组更加准确的变换关系(六常数模型)。为从不变性元组中找出相似三角形(五角星和圆圈距离最近的点)，我们使用k-d树进行最近邻搜索，搜索的时间复杂度为O(nlogm)，其中，n为列表I生成的不变性元组数量，m为列表R生成的不变性元组数量。具体细节可参考Astroalign(https://pypi.org/project/astroalign/)图像对齐算法。为适应本文星像匹配算法，我们对Astroalign默认设置的参数进行微调。

1.4 使用k-means聚类算法求局部变换

使用以上方法即可求出一个正确的变换，从而将列表I中星像的像素坐标转换为天球坐标。但是，当图像较大时，由于图像扭曲效应，一个变换并不能满足整张图像中星像像素坐标向天球坐标的转换。由于图像扭曲的局部效应差别不大，本文考虑使用k-means聚类算法。其基本思想是以空间中的k个点为中心，对距离k个中心点最近的对象分别归类。通过迭代方法，逐次更新各聚类中心的值，直到得到满足要求的聚类结果。该算法可以将距离较近的星像归为一类，然后对每一类星像单独求变换，简称局部变换。在整体变换的基础上，我们很容易找到图像中每一类星像在列表R中相对应的子区域，进而重新求解一个变换作为局部变换。当某一类星像由于数量过少或亮星数量不够时，可能无法求得局部变换，此时使用整体变换代替局部变换。

1.5 最终验证

我们使用整体或局部变换将列表I中所有星像的像素坐标转换为天球坐标，在列表R中找出和列表I中星像最近的星像。使用k-d树的最近邻搜索算法和文[13]中的距离公式为

s2=(Δαcosδ)2+(Δδ)2，

(3)

其中，s为距离的偏差；Δα为赤经方向的偏差；Δδ为赤纬方向的偏差；δ为列表I中星像的赤纬。使用平方公式可以减少开根号的时间开销。由于大气折射较差的影响可以达到约2″(在天顶距为60°,视场为1°的极端情况下)，再加上视场畸变等因素，通常，当s小于3″的阈值时(可调)，认为列表I和列表R中的星像能够相匹配。选用较大的阈值可以匹配更多的星像，以便做进一步分析。

2 实验分析

为验证本文算法的可行性，我们使用云南天文台1 m光学望远镜拍摄的稀疏星场图像(星像数量～100颗，图像大小为4 096×411 2，视场中心的赤经为～37.162°，赤纬为～-8.339°)和2.4 m光学望远镜拍摄的密集星场图像(星像数量～1 000颗，图像大小为1 900×1 900，视场中心的赤经为～72.791°，赤纬为～+43.680°)进行实验。图3展示了其中两幅图像的匹配情况。其中，1 m光学望远镜的比例尺为～0.″234/pixel，2.4m光学望远镜的比例尺为～0.″286/pixel。本文方法是对文[8]基于三角形方法的一种改进。由于不同的人所使用的机器不同或者使用的编程语言不一样，加上亮星的数量不同等原因，所以绝对的运行时间不具有可比性。

图3 (a)使用云南天文台1 m光学望远镜拍摄的稀疏星场图像匹配情况；(b)使用云南天文台2.4 m光学望远镜拍摄的密集星场图像匹配情况(红色圆圈内的星像表示搜星成功，但是没有匹配上的星像。为尽可能地将暗星匹配上，搜星程序选用较低的检测阈值，所以出现较多的红色圆圈内的伪检测对象)Fig.3 (a)Matching of sparse star field images captured using the 1 m optical telescope at the Yunnan Observatories；(b)matching of dense star field images captured using the 2.4 m optical telescope at the Yunnan Observatories (The stars in the red circle indicate successful search，but there are no matching stars. To match faint stars as much as possible，the star search program uses a lower detection threshold，resulting in more false detection objects within the red circle)

对于其中一幅稀疏星场图像，本文算法求得的图像比例尺为0.″234 10/pixel，搜星数量为71颗，最终匹配成功62颗，匹配率为87.324%。使用(3)式开根号计算残差，求得平均残差为0.″121。若使用k-means聚类算法对图像中的星像进行分类，例如分为4类，分别求图像的局部变换。使用局部变换来做星像匹配，求得图像的比例尺分别为0.″234 04/pixel，0.″234 10/pixel，0.″234 00/pixel和0.″234 19/pixel。星像匹配数量不变，但是平均残差为0.″105，平均残差明显降低。

为进一步验证使用k-means聚类算法的优势，我们使用云南天文台2.4 m光学望远镜拍摄的密集星场图像进行实验，搜星数量为1 290。当使用整体变换时，匹配数量934颗，匹配率为72.403%，平均残差为0.″202。当使用局部变换时，例如使用k-means聚类算法将整幅图像分为1，4，9，…，100个小区域时，平均残差和局部变换匹配率(局部变换数量和分类数量的比值)如图4。当分为49类时，所得平均残差最小，为0.″173。当分类继续增多时，平均残差不再减小，我们猜测主要原因是找不到足够的局部变换，只能使用整体变换代替。当分为64类时，找到的局部变换数量增多不再明显，局部变换匹配率降低。因此，并不是分类越多，平均残差就越小。在此次实验中，匹配上的星像数量变化不大：分为1和81类的情况下，匹配星像数量为934颗；分为4，9，16，36，49和100类的情况下，匹配星像数量为935颗；分为25类的情况下，匹配星像数量为936颗；分为64类的情况下，匹配星像数量为928颗。此外，在求局部变换时，所求得的比例尺也有微小差别：当分区较少时，比如4个分区时，比例尺差异较小，在万分之几角秒每像素；当分区较多时，比如49个分区时，比例尺差异在万分之几到千分之几角秒每像素。总体来说，分区距离越近，比例尺的差异越小；分区距离越远，比例尺的差异越大。图5展示了该密集星场图像使用k-means聚类算法分为49类时的图像比例尺情况。

图4 不同分类数量的平均残差与局部变换匹配率Fig.4 Average residual and local transformation matching rate of different number of clusters

注：单位为角秒/像素；“None”表示求局部变换失败。图5 使用k-means聚类算法分为49类时的图像比例尺情况Fig.5 Image scale when using k-means clustering algorithm to divide into 49 categories

实验表明，本文提出的星像匹配算法是一种可行的方案，且本文使用k-means聚类算法对图像进行分割，在一定程度上可以提高匹配的精度，同时能够自适应图像局部比例尺的变化。

3 解释与讨论

本文对文[8]的星像匹配方法进行了改进。文[8]的方法依次遍历搜索最近邻点，其时间复杂度为O(nm)。例如，对于列表I中星像的数量为1 290，列表R中星像的数量为12 634的情况，从列表R中找出距离列表I中星像最近的星像，使用k-d树实现，只需要1.288 s。如果使用依次遍历搜索算法(Python实现)，耗时超过50 s。k-d树是一种空间划分树，适用于最近邻搜索，其时间复杂度为O(nlogm)。

对于求局部变换来说，有一种规则分区的方法是将图像分割成M×N个矩形区域，虽然速度快，但是规则分区方法不能很好地适应星像密度分布不均匀的情况，可能在某些分区内有过多的星像，在另一些分区内星像数目太少，甚至没有。mean-shift和DBSCAN(Density-Based Spatial Clustering of Applications with Noise)虽然也是常用的聚类算法，但事实上，使用不同的聚类方法并不影响匹配的精度。基于本文的实验资料，k-means算法的聚类效果具有良好的适用性。

在k-means聚类算法分区的实际应用中，需要按k=M×N的形式遍历k的值，以确保每个分区类似正方形或圆形，从而解决视场畸变情况下的配准问题。理论上，每个分区至少需要3颗星(以便求解四常数模型)才能进行配准。但从代码实现上，需要至少5颗星(生成不变性元组时，以5颗相距最近的星为一组)才能进行配准。

本文使用的RANSAC算法来源于Astroalign。该算法有两个参数阈值需要调整：第1个是可容忍的全等三角形匹配像素误差；第2个是需要匹配的全等三角形数量。设置第1个参数是由于误差的影响，两个三角形无论怎么变换，也很难完全相同。设置一个可容忍的全等三角形匹配像素误差阈值是为了匹配全等三角形，该阈值是经验性的，通常为2个像素。设置过小，全等三角形无法匹配，设置过大，会出现误匹配。实际应用还需要根据图像的分辨率大小确定。调整第2个参数阈值是为了让匹配的相似三角形经过变换后，如果满足全等匹配的三角形的数量超过这个阈值，则认为找到了一个正确的变换。理论上，该参数设置越大越好，但是由于三角形全等匹配存在一定误差，需要设置稍微小一点。经实验验证，在找不到一个正确变换时，可以尝试将此阈值设置小一点，以此来换取一个正确的变换。然而，当此阈值太小时，比如为3时，可能出现找到一个错误的正确变换而无法应用于最终匹配阶段的情况。因此，该阈值应从大到小设置，典型地，可以从10开始递减设置。

4 结 论

本文提出了一种基于k-d树和k-means聚类算法的星像匹配算法，该算法适用于具有明显畸变的大视场。本文方法的优势主要体现在以下两方面：(1)本文灵活地调用Python库函数，更方便实现；(2)本文匹配算法和空间存储的优化。相对于文[8]的方法，本文采用了Astroalign中的三角形不变性元组进行相似三角形的匹配，不需要输入比例尺参数。而文[8]的方法采用全等三角形的方式进行匹配，需要预先输入望远镜精确的比例尺参数。此外，使用k-d树进行最近邻搜索能够有效地减少时间复杂度，使用k-means聚类能够有效避免图像中部分区域没有星像的情况。

致谢：感谢云南天文台1 m望远镜和2.4 m望远镜工作人员的支持。