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

(1. 暨南大學(xué)計(jì)算機(jī)科學(xué)系，廣東 廣州 510632；2. 暨南大學(xué)中法天體測(cè)量、動(dòng)力學(xué)與空間科學(xué)聯(lián)合實(shí)驗(yàn)室，廣東 廣州 510632；3. 暨南大學(xué)物理學(xué)系，廣東 廣州 510632)

在天文觀測(cè)中，為了方便資料的歸算和處理，通常需要將圖像中星像的位置和光度等信息與參考列表中對(duì)應(yīng)的信息進(jìn)行匹配，簡(jiǎn)稱星像匹配。常見(jiàn)的星像匹配方法有基于向量的方法[1]，基于徑向和環(huán)向特征的方法[2]，基于三角形的方法[3]和基于局部敏感哈希的方法[4]等。比較成功的應(yīng)用是Astrometry.net[5]，它使用四邊形代替三角形，并通過(guò)貝葉斯決策準(zhǔn)則來(lái)選擇正確的圖像變換。文[6]研究發(fā)現(xiàn)，基于向量的方法比基于徑向和環(huán)向特征的方法更優(yōu)，具體表現(xiàn)為基于向量的方法耗時(shí)更少，匹配率更高。雖然一些研究指出，四邊形匹配算法比三角形匹配算法具有更好的效果[4]。但是具有簡(jiǎn)單結(jié)構(gòu)的三角形仍是較為常用的一種星像匹配算法[1，3，7]。文[8]注意到，基于三角形的方法是以上3種方法[1-3]中匹配率最高，即最穩(wěn)定的一種方法。但是，該方法列表中的星所構(gòu)成的三角形數(shù)目巨大，造成存儲(chǔ)和檢索的困難。不僅如此，匹配需要構(gòu)造全等三角形，匹配全等三角形需要預(yù)先知道圖像的比例尺。圖像的比例尺不是一成不變的，受到各種因素的影響，如焦距、溫度、圖像扭曲和鏡筒彎曲等。

受Astroalign[9]圖像對(duì)齊算法的啟發(fā)，本文提出一種基于k-d樹(shù)和k-means聚類(lèi)算法的星像匹配算法。該算法利用k-d樹(shù)進(jìn)行三角形盲匹配以找到整幅圖像區(qū)域的整體變換，然后通過(guò)k-means聚類(lèi)算法將整幅圖像分割為多個(gè)小區(qū)域，從而求解小區(qū)域范圍內(nèi)的局部變換。如果不存在這樣的局部變換，則用整體變換代替。其中，使用k-d樹(shù)的最近鄰搜索算法可以提高匹配的速度，使用k-means聚類(lèi)算法可以提高匹配的精度。該算法可以間接計(jì)算圖像的比例尺，從而自適應(yīng)地體現(xiàn)圖像比例尺的微小變化。

1 方 法

因?yàn)榍叭说姆椒ㄊ腔谌热切蔚钠ヅ鋄1]，所以需要預(yù)先知道圖像的比例尺。雖然獲得圖像的比例尺并不困難，但是圖像的比例尺會(huì)因各種因素的影響而發(fā)生微小的變化。因此，理論上，星像匹配應(yīng)使用相似三角形而非全等三角形。為了減小圖像比例尺變化和生成大量三角形的影響，本文3次使用k-d樹(shù)的最近鄰搜索算法：第1次使用k-d樹(shù)生成三角形不變性元組；第2次使用k-d樹(shù)匹配三角形不變性元組，從而匹配相似三角形；第3次使用k-d樹(shù)進(jìn)行最終驗(yàn)證。本質(zhì)上，k-d樹(shù)是一種平衡二叉樹(shù)，也是一種空間劃分樹(shù)。它使用二分的方式將整個(gè)k維向量空間不斷劃分為若干個(gè)局部空間。該算法能夠在每次搜索過(guò)程中進(jìn)行分支判斷，選擇滿足要求的局部子空間。該方法有效避免了全局空間搜索，以達(dá)到快速匹配的目的。為提高匹配的精度，我們使用k-means聚類(lèi)算法進(jìn)行區(qū)域劃分，以達(dá)到局部精準(zhǔn)匹配的目的，使用心射投影(中心投影)進(jìn)行天球坐標(biāo)系到標(biāo)準(zhǔn)坐標(biāo)系的投影。在不考慮投影誤差的前提下，本文假定參考列表的天球坐標(biāo)投影到了標(biāo)準(zhǔn)坐標(biāo)，同時(shí)也知道了圖像中星像的像素坐標(biāo)(參考星的數(shù)目與圖像中星像的數(shù)目不一定相等)。我們的主要任務(wù)是找出共同對(duì)象中參考星的標(biāo)準(zhǔn)坐標(biāo)和圖像中星像像素坐標(biāo)的一一對(duì)應(yīng)關(guān)系。

1.1 創(chuàng)建星像列表和讀取參考列表

創(chuàng)建星像列表需要對(duì)天文圖像進(jìn)行搜星和測(cè)光。本文使用DAOFIND[10]程序進(jìn)行搜星與測(cè)光，然后根據(jù)星像的亮度進(jìn)行降序排序，獲得星像列表I。列表I中至少包含像素坐標(biāo)和光通量?jī)煞N信息。根據(jù)天文圖像頭文件中的視場(chǎng)中心、圖像大小和焦距等信息推斷圖像對(duì)應(yīng)的天球坐標(biāo)的視場(chǎng)。從Gaia DR3[11]中下載對(duì)應(yīng)的視場(chǎng)中的星像，并按照星等從小到大排序獲得參考列表R，列表R的視場(chǎng)遠(yuǎn)遠(yuǎn)大于列表I。參考列表R至少包含天球坐標(biāo)(赤經(jīng)和赤緯)、自行(實(shí)際計(jì)算中可以依據(jù)觀測(cè)歷元和自行計(jì)算出星表在觀測(cè)歷元的近似位置)和星等3種信息。

1.2 生成不變性元組

本文使用相似三角形進(jìn)行三角形盲匹配。一對(duì)匹配好的相似三角形可以求解一個(gè)變換。由于匹配的是相似三角形，所以本文算法還可以間接地求解CCD圖像的比例尺。為對(duì)相似三角形進(jìn)行盲匹配，本文使用Astroalign算法中的三角形不變性元組。假設(shè)三角形的三條邊分別為L(zhǎng)0，L1和L2，三角形不變性元組可以表示為

(1)

(2)

為了分析三角形不變性元組的特性，直觀地了解三角形不變性元組，我們使用天文圖像和參考列表中的星像構(gòu)建三角形來(lái)生成不變性元組，并繪制如圖1的三角形不變性元組特征圖。生成三角形不變性元組的流程如圖2。其中，兩次使用了k-d樹(shù)，第1次使用k-d樹(shù)的目的是找出5個(gè)在位置上最近的點(diǎn)，第2次使用k-d樹(shù)的目的是通過(guò)不變性元組找出可能匹配的三角形。選用5個(gè)位置最近的點(diǎn)為一組的分組，可以減少三角形生成的數(shù)量，且有利于局部匹配。因?yàn)闄M跨小片區(qū)域的三角形和橫跨大片區(qū)域的三角形即使相似，也不可能導(dǎo)出一個(gè)正確的變換，而且橫跨大片區(qū)域的三角形還可能存在圖像扭曲。第2次使用k-d樹(shù)可以快速找出匹配的相似三角形，從而嘗試求解一個(gè)變換。

圖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 導(dǎo)出變換

由圖1可知，五角星(☆)和圓圈(○)距離越近，表示兩個(gè)三角形相似度越高。兩個(gè)相似三角形的頂點(diǎn)對(duì)應(yīng)關(guān)系可以嘗試求解一組變換關(guān)系。本文使用Astroalign中的RANSAC[12]算法，主要體現(xiàn)為如果星像列表I中的星像進(jìn)行該變換后，部分三角形能夠滿足全等匹配關(guān)系，說(shuō)明該變換是一個(gè)星像列表I轉(zhuǎn)換到參考列表R的正確變換；否則，重新選擇兩個(gè)相似三角形求解另一組變換關(guān)系。重復(fù)以上步驟，直到獲得一組正確的變換關(guān)系。當(dāng)找到一組正確的變換關(guān)系后，即可找出亮星在像素坐標(biāo)和標(biāo)準(zhǔn)坐標(biāo)之間的對(duì)應(yīng)關(guān)系，最后使用最小二乘擬合出一組更加準(zhǔn)確的變換關(guān)系(六常數(shù)模型)。為從不變性元組中找出相似三角形(五角星和圓圈距離最近的點(diǎn))，我們使用k-d樹(shù)進(jìn)行最近鄰搜索，搜索的時(shí)間復(fù)雜度為O(nlogm)，其中，n為列表I生成的不變性元組數(shù)量，m為列表R生成的不變性元組數(shù)量。具體細(xì)節(jié)可參考Astroalign(https://pypi.org/project/astroalign/)圖像對(duì)齊算法。為適應(yīng)本文星像匹配算法，我們對(duì)Astroalign默認(rèn)設(shè)置的參數(shù)進(jìn)行微調(diào)。

1.4 使用k-means聚類(lèi)算法求局部變換

使用以上方法即可求出一個(gè)正確的變換，從而將列表I中星像的像素坐標(biāo)轉(zhuǎn)換為天球坐標(biāo)。但是，當(dāng)圖像較大時(shí)，由于圖像扭曲效應(yīng)，一個(gè)變換并不能滿足整張圖像中星像像素坐標(biāo)向天球坐標(biāo)的轉(zhuǎn)換。由于圖像扭曲的局部效應(yīng)差別不大，本文考慮使用k-means聚類(lèi)算法。其基本思想是以空間中的k個(gè)點(diǎn)為中心，對(duì)距離k個(gè)中心點(diǎn)最近的對(duì)象分別歸類(lèi)。通過(guò)迭代方法，逐次更新各聚類(lèi)中心的值，直到得到滿足要求的聚類(lèi)結(jié)果。該算法可以將距離較近的星像歸為一類(lèi)，然后對(duì)每一類(lèi)星像單獨(dú)求變換，簡(jiǎn)稱局部變換。在整體變換的基礎(chǔ)上，我們很容易找到圖像中每一類(lèi)星像在列表R中相對(duì)應(yīng)的子區(qū)域，進(jìn)而重新求解一個(gè)變換作為局部變換。當(dāng)某一類(lèi)星像由于數(shù)量過(guò)少或亮星數(shù)量不夠時(shí)，可能無(wú)法求得局部變換，此時(shí)使用整體變換代替局部變換。

1.5 最終驗(yàn)證

我們使用整體或局部變換將列表I中所有星像的像素坐標(biāo)轉(zhuǎn)換為天球坐標(biāo)，在列表R中找出和列表I中星像最近的星像。使用k-d樹(shù)的最近鄰搜索算法和文[13]中的距離公式為

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

(3)

其中，s為距離的偏差；Δα為赤經(jīng)方向的偏差；Δδ為赤緯方向的偏差；δ為列表I中星像的赤緯。使用平方公式可以減少開(kāi)根號(hào)的時(shí)間開(kāi)銷(xiāo)。由于大氣折射較差的影響可以達(dá)到約2″(在天頂距為60°,視場(chǎng)為1°的極端情況下)，再加上視場(chǎng)畸變等因素，通常，當(dāng)s小于3″的閾值時(shí)(可調(diào))，認(rèn)為列表I和列表R中的星像能夠相匹配。選用較大的閾值可以匹配更多的星像，以便做進(jìn)一步分析。

2 實(shí)驗(yàn)分析

為驗(yàn)證本文算法的可行性，我們使用云南天文臺(tái)1 m光學(xué)望遠(yuǎn)鏡拍攝的稀疏星場(chǎng)圖像(星像數(shù)量～100顆，圖像大小為4 096×411 2，視場(chǎng)中心的赤經(jīng)為～37.162°，赤緯為～-8.339°)和2.4 m光學(xué)望遠(yuǎn)鏡拍攝的密集星場(chǎng)圖像(星像數(shù)量～1 000顆，圖像大小為1 900×1 900，視場(chǎng)中心的赤經(jīng)為～72.791°，赤緯為～+43.680°)進(jìn)行實(shí)驗(yàn)。圖3展示了其中兩幅圖像的匹配情況。其中，1 m光學(xué)望遠(yuǎn)鏡的比例尺為～0.″234/pixel，2.4m光學(xué)望遠(yuǎn)鏡的比例尺為～0.″286/pixel。本文方法是對(duì)文[8]基于三角形方法的一種改進(jìn)。由于不同的人所使用的機(jī)器不同或者使用的編程語(yǔ)言不一樣，加上亮星的數(shù)量不同等原因，所以絕對(duì)的運(yùn)行時(shí)間不具有可比性。

圖3 (a)使用云南天文臺(tái)1 m光學(xué)望遠(yuǎn)鏡拍攝的稀疏星場(chǎng)圖像匹配情況；(b)使用云南天文臺(tái)2.4 m光學(xué)望遠(yuǎn)鏡拍攝的密集星場(chǎng)圖像匹配情況(紅色圓圈內(nèi)的星像表示搜星成功，但是沒(méi)有匹配上的星像。為盡可能地將暗星匹配上，搜星程序選用較低的檢測(cè)閾值，所以出現(xiàn)較多的紅色圓圈內(nèi)的偽檢測(cè)對(duì)象)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)

對(duì)于其中一幅稀疏星場(chǎng)圖像，本文算法求得的圖像比例尺為0.″234 10/pixel，搜星數(shù)量為71顆，最終匹配成功62顆，匹配率為87.324%。使用(3)式開(kāi)根號(hào)計(jì)算殘差，求得平均殘差為0.″121。若使用k-means聚類(lèi)算法對(duì)圖像中的星像進(jìn)行分類(lèi)，例如分為4類(lèi)，分別求圖像的局部變換。使用局部變換來(lái)做星像匹配，求得圖像的比例尺分別為0.″234 04/pixel，0.″234 10/pixel，0.″234 00/pixel和0.″234 19/pixel。星像匹配數(shù)量不變，但是平均殘差為0.″105，平均殘差明顯降低。

為進(jìn)一步驗(yàn)證使用k-means聚類(lèi)算法的優(yōu)勢(shì)，我們使用云南天文臺(tái)2.4 m光學(xué)望遠(yuǎn)鏡拍攝的密集星場(chǎng)圖像進(jìn)行實(shí)驗(yàn)，搜星數(shù)量為1 290。當(dāng)使用整體變換時(shí)，匹配數(shù)量934顆，匹配率為72.403%，平均殘差為0.″202。當(dāng)使用局部變換時(shí)，例如使用k-means聚類(lèi)算法將整幅圖像分為1，4，9，…，100個(gè)小區(qū)域時(shí)，平均殘差和局部變換匹配率(局部變換數(shù)量和分類(lèi)數(shù)量的比值)如圖4。當(dāng)分為49類(lèi)時(shí)，所得平均殘差最小，為0.″173。當(dāng)分類(lèi)繼續(xù)增多時(shí)，平均殘差不再減小，我們猜測(cè)主要原因是找不到足夠的局部變換，只能使用整體變換代替。當(dāng)分為64類(lèi)時(shí)，找到的局部變換數(shù)量增多不再明顯，局部變換匹配率降低。因此，并不是分類(lèi)越多，平均殘差就越小。在此次實(shí)驗(yàn)中，匹配上的星像數(shù)量變化不大：分為1和81類(lèi)的情況下，匹配星像數(shù)量為934顆；分為4，9，16，36，49和100類(lèi)的情況下，匹配星像數(shù)量為935顆；分為25類(lèi)的情況下，匹配星像數(shù)量為936顆；分為64類(lèi)的情況下，匹配星像數(shù)量為928顆。此外，在求局部變換時(shí)，所求得的比例尺也有微小差別：當(dāng)分區(qū)較少時(shí)，比如4個(gè)分區(qū)時(shí)，比例尺差異較小，在萬(wàn)分之幾角秒每像素；當(dāng)分區(qū)較多時(shí)，比如49個(gè)分區(qū)時(shí)，比例尺差異在萬(wàn)分之幾到千分之幾角秒每像素?？傮w來(lái)說(shuō)，分區(qū)距離越近，比例尺的差異越?。环謪^(qū)距離越遠(yuǎn)，比例尺的差異越大。圖5展示了該密集星場(chǎng)圖像使用k-means聚類(lèi)算法分為49類(lèi)時(shí)的圖像比例尺情況。

圖4 不同分類(lèi)數(shù)量的平均殘差與局部變換匹配率Fig.4 Average residual and local transformation matching rate of different number of clusters

注：?jiǎn)挝粸榻敲?像素；“None”表示求局部變換失敗。圖5 使用k-means聚類(lèi)算法分為49類(lèi)時(shí)的圖像比例尺情況Fig.5 Image scale when using k-means clustering algorithm to divide into 49 categories

實(shí)驗(yàn)表明，本文提出的星像匹配算法是一種可行的方案，且本文使用k-means聚類(lèi)算法對(duì)圖像進(jìn)行分割，在一定程度上可以提高匹配的精度，同時(shí)能夠自適應(yīng)圖像局部比例尺的變化。

3 解釋與討論

本文對(duì)文[8]的星像匹配方法進(jìn)行了改進(jìn)。文[8]的方法依次遍歷搜索最近鄰點(diǎn)，其時(shí)間復(fù)雜度為O(nm)。例如，對(duì)于列表I中星像的數(shù)量為1 290，列表R中星像的數(shù)量為12 634的情況，從列表R中找出距離列表I中星像最近的星像，使用k-d樹(shù)實(shí)現(xiàn)，只需要1.288 s。如果使用依次遍歷搜索算法(Python實(shí)現(xiàn))，耗時(shí)超過(guò)50 s。k-d樹(shù)是一種空間劃分樹(shù)，適用于最近鄰搜索，其時(shí)間復(fù)雜度為O(nlogm)。

對(duì)于求局部變換來(lái)說(shuō)，有一種規(guī)則分區(qū)的方法是將圖像分割成M×N個(gè)矩形區(qū)域，雖然速度快，但是規(guī)則分區(qū)方法不能很好地適應(yīng)星像密度分布不均勻的情況，可能在某些分區(qū)內(nèi)有過(guò)多的星像，在另一些分區(qū)內(nèi)星像數(shù)目太少，甚至沒(méi)有。mean-shift和DBSCAN(Density-Based Spatial Clustering of Applications with Noise)雖然也是常用的聚類(lèi)算法，但事實(shí)上，使用不同的聚類(lèi)方法并不影響匹配的精度。基于本文的實(shí)驗(yàn)資料，k-means算法的聚類(lèi)效果具有良好的適用性。

在k-means聚類(lèi)算法分區(qū)的實(shí)際應(yīng)用中，需要按k=M×N的形式遍歷k的值，以確保每個(gè)分區(qū)類(lèi)似正方形或圓形，從而解決視場(chǎng)畸變情況下的配準(zhǔn)問(wèn)題。理論上，每個(gè)分區(qū)至少需要3顆星(以便求解四常數(shù)模型)才能進(jìn)行配準(zhǔn)。但從代碼實(shí)現(xiàn)上，需要至少5顆星(生成不變性元組時(shí)，以5顆相距最近的星為一組)才能進(jìn)行配準(zhǔn)。

本文使用的RANSAC算法來(lái)源于Astroalign。該算法有兩個(gè)參數(shù)閾值需要調(diào)整：第1個(gè)是可容忍的全等三角形匹配像素誤差；第2個(gè)是需要匹配的全等三角形數(shù)量。設(shè)置第1個(gè)參數(shù)是由于誤差的影響，兩個(gè)三角形無(wú)論怎么變換，也很難完全相同。設(shè)置一個(gè)可容忍的全等三角形匹配像素誤差閾值是為了匹配全等三角形，該閾值是經(jīng)驗(yàn)性的，通常為2個(gè)像素。設(shè)置過(guò)小，全等三角形無(wú)法匹配，設(shè)置過(guò)大，會(huì)出現(xiàn)誤匹配。實(shí)際應(yīng)用還需要根據(jù)圖像的分辨率大小確定。調(diào)整第2個(gè)參數(shù)閾值是為了讓匹配的相似三角形經(jīng)過(guò)變換后，如果滿足全等匹配的三角形的數(shù)量超過(guò)這個(gè)閾值，則認(rèn)為找到了一個(gè)正確的變換。理論上，該參數(shù)設(shè)置越大越好，但是由于三角形全等匹配存在一定誤差，需要設(shè)置稍微小一點(diǎn)。經(jīng)實(shí)驗(yàn)驗(yàn)證，在找不到一個(gè)正確變換時(shí)，可以嘗試將此閾值設(shè)置小一點(diǎn)，以此來(lái)?yè)Q取一個(gè)正確的變換。然而，當(dāng)此閾值太小時(shí)，比如為3時(shí)，可能出現(xiàn)找到一個(gè)錯(cuò)誤的正確變換而無(wú)法應(yīng)用于最終匹配階段的情況。因此，該閾值應(yīng)從大到小設(shè)置，典型地，可以從10開(kāi)始遞減設(shè)置。

4 結(jié) 論

本文提出了一種基于k-d樹(shù)和k-means聚類(lèi)算法的星像匹配算法，該算法適用于具有明顯畸變的大視場(chǎng)。本文方法的優(yōu)勢(shì)主要體現(xiàn)在以下兩方面：(1)本文靈活地調(diào)用Python庫(kù)函數(shù)，更方便實(shí)現(xiàn)；(2)本文匹配算法和空間存儲(chǔ)的優(yōu)化。相對(duì)于文[8]的方法，本文采用了Astroalign中的三角形不變性元組進(jìn)行相似三角形的匹配，不需要輸入比例尺參數(shù)。而文[8]的方法采用全等三角形的方式進(jìn)行匹配，需要預(yù)先輸入望遠(yuǎn)鏡精確的比例尺參數(shù)。此外，使用k-d樹(shù)進(jìn)行最近鄰搜索能夠有效地減少時(shí)間復(fù)雜度，使用k-means聚類(lèi)能夠有效避免圖像中部分區(qū)域沒(méi)有星像的情況。

致謝：感謝云南天文臺(tái)1 m望遠(yuǎn)鏡和2.4 m望遠(yuǎn)鏡工作人員的支持。