基于共有星的快速星圖識(shí)別方法

楊 上，周召發(fā)，劉先一

（火箭軍工程大學(xué)兵器發(fā)射理論與技術(shù)國(guó)家重點(diǎn)學(xué)科實(shí)驗(yàn)室，陜西西安710025）

1 引 言

數(shù)字天頂儀是一種高精度天文定位儀器，在大地天文測(cè)量領(lǐng)域得到廣泛應(yīng)用［1-4］。星圖識(shí)別為數(shù)字天頂儀定位提供起算數(shù)據(jù)，傳統(tǒng)數(shù)字天頂儀的星圖識(shí)別是采用三角形星圖識(shí)別與坐標(biāo)轉(zhuǎn)換結(jié)合方法對(duì)旋轉(zhuǎn)拍攝的16幅星圖進(jìn)行識(shí)別［5］。經(jīng)大量實(shí)驗(yàn)發(fā)現(xiàn)，數(shù)字天頂儀在一次定位循環(huán)拍攝的16幅星圖中存在著大量的共有星，即同一恒星會(huì)成像在多幅星圖中的情況，若對(duì)拍攝的16幅星圖上所有星點(diǎn)都單獨(dú)進(jìn)行識(shí)別，增加了星圖識(shí)別的冗余量，降低了星圖識(shí)別的效率，進(jìn)而增加了數(shù)字天頂儀的定位所需的時(shí)間。對(duì)此本文提出了一種基于共有星的快速星圖識(shí)別方法。該方法的核心思想首先對(duì)一個(gè)定位循環(huán)中在位置9和位置10拍攝的星圖單獨(dú)進(jìn)行三角形星圖識(shí)別，其次建立星點(diǎn)歸算模型，找出剩余位置的星圖與已識(shí)別星圖間的共有星，最后對(duì)剩余星圖上新增星單獨(dú)進(jìn)行識(shí)別，以此提高星圖識(shí)別的速度。

本文首先介紹了單幅星圖的三角形識(shí)別原理，其次通過構(gòu)造輔助坐標(biāo)系建立了傾斜修正模型，消除精調(diào)平狀態(tài)下的小傾角，再在此基礎(chǔ)上建立了星點(diǎn)歸算模型，找出星圖間的共有星，最后通過實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證該星圖識(shí)別方法的合理性。

2 數(shù)字天頂儀星圖識(shí)別原理

數(shù)字天頂儀進(jìn)行星圖識(shí)別的導(dǎo)航星表是依巴谷星表，依據(jù)導(dǎo)航星表可得到恒星的赤經(jīng)赤緯信息，再經(jīng)視位置計(jì)算和時(shí)間補(bǔ)償可得到恒星的天文經(jīng)度α、天文緯度δ，則恒星的天球切平面坐標(biāo)（ξ，η）可按下式計(jì)算：

其中，（α0，δ0）為GPS提供的概略天文經(jīng)緯度。

恒星的理論CCD像點(diǎn)坐標(biāo)為：

恒星的理論像距可按下式計(jì)算：

理想情況下，CCD平面為水平面，恒星星光通過數(shù)字天頂儀的焦距成像在CCD平面上，假設(shè)像點(diǎn)坐標(biāo)為 （xm，ym） ，（xn，yn） ，則恒星像點(diǎn)像距按下式計(jì)算：

若恒星理論像點(diǎn)組成的三角形與恒星像點(diǎn)組成的三角形滿足下式則匹配成功：

其中，ε為測(cè)量不確定度。

若對(duì)數(shù)字天頂儀拍攝的星圖都采用三角形星圖識(shí)別，則計(jì)算量較大，原因在于觀測(cè)三角形與導(dǎo)航三角形容量較大［6］，降低了數(shù)字天頂儀星圖識(shí)別效率，對(duì)此數(shù)字天頂儀采用三角形識(shí)別與坐標(biāo)轉(zhuǎn)換相結(jié)合的方式進(jìn)行星圖識(shí)別，識(shí)別過程如下：

（1）采用三角形星圖識(shí)別原理識(shí)別出6顆亮星。

（2）利用已識(shí)別出亮星的理論像點(diǎn)坐標(biāo)與CCD像點(diǎn)坐標(biāo)，建立仿射變換模型［7］，如下式：

其中，仿射變換模型的參數(shù)可由最小二乘法計(jì)算得到。

（3）計(jì)算數(shù)字天頂儀視場(chǎng)內(nèi)恒星的理論像點(diǎn)坐標(biāo)。

（4）利用式（6）計(jì)算出轉(zhuǎn)換后的恒星像點(diǎn)坐標(biāo)（x′i，yi′） 。

（5）將轉(zhuǎn)換后的像點(diǎn)坐標(biāo)與CCD像點(diǎn)坐標(biāo)進(jìn)行匹配，若滿足下式則匹配成功：

3 傾斜修正

傳統(tǒng)的數(shù)字天頂儀星圖識(shí)別是在精密整平狀態(tài)下進(jìn)行的，認(rèn)為CCD平面是水平面，但即使在精調(diào)平狀態(tài)下其CCD平面與水平面仍存在±10″內(nèi)的小傾角［8］，由下文第4節(jié)可知，建立的星點(diǎn)歸算模型與CCD北向基準(zhǔn)坐標(biāo)有關(guān)，而北向基準(zhǔn)坐標(biāo)系的建立與水平面為基準(zhǔn)［2］，為此需要修正實(shí)際拍攝狀態(tài)下存在的小傾角。

實(shí)際拍攝狀態(tài)下，由文獻(xiàn)［9］和文獻(xiàn)［10］可知，經(jīng)過較標(biāo)的傾角儀敏感軸組成的平面與CCD平面平行，故可通過傾角儀測(cè)出CCD圖像坐標(biāo)系兩坐標(biāo)軸與水平面之間的夾角，假設(shè)CCD圖像坐標(biāo)系兩坐標(biāo)軸與水平狀面之間的夾角θ1，θ2，可按下式計(jì)算得到：

式中，m，n為傾角儀的輸出值；φ為傾角儀m軸與CCD圖像坐標(biāo)系x軸之間的夾角。

如圖1所示，實(shí)線框?yàn)榕臄z位置的水平面，坐標(biāo)系O-xTyTzT為實(shí)際拍攝狀態(tài)下CCD圖像坐標(biāo)系，其中CCD圖像坐標(biāo)系的坐標(biāo)原點(diǎn)在CCD敏感器的中心，兩坐標(biāo)軸分別與CCD敏感器的兩條邊平行。OC，OD為OxT軸與OyT軸在水平面上的投影，OC，OD不一定垂直，則∠COxT=θ1，∠DOyT=θ2，AB為實(shí)際拍攝狀態(tài)下的CCD平面與水平面之間的交線。

圖1 CCD傾斜修正Fig.1 CCD tilt correct

由二面角的相關(guān)定理，可得實(shí)際拍攝狀態(tài)下的CCD平面與水平面之間構(gòu)成的二面角θ為：

如圖1所示，構(gòu)造輔助坐標(biāo)系O-xFyFzF，其中OxF軸在實(shí)際拍攝狀態(tài)下的CCD平面內(nèi)且OxF⊥AB，OyF軸沿AB方向，OzF軸與OxF軸和OyF軸符合右手定則，OxS軸為 OxF軸的投影，則∠xFOxS=θ。設(shè)∠xFOxT=α，由二面角的相關(guān)定理可得：

根據(jù)方向余弦坐標(biāo)變換原理可得存在方向余弦矩陣C1，可將坐標(biāo)系O-xTyTzT中的星點(diǎn)坐標(biāo)轉(zhuǎn)換到坐標(biāo)系O-xFyFzF中，即：

以水平面為基準(zhǔn)，建立水平狀態(tài)下的坐標(biāo)系O-xSySzS，其中OxS軸為OxF軸的投影，OyS軸與OyF軸重合，OzS軸與OxS軸和OyS軸滿足右手定則，則輔助坐標(biāo)系O-xFyFzF可繞OyF旋轉(zhuǎn)θ角與坐標(biāo)系O-xSySzS重合，則在輔助坐標(biāo)系O-xFyFzF中，存在方向矢量 （xF，yF， -f） 滿足下式：

恒星位于無窮遠(yuǎn)處，發(fā)出的光為平行光，假設(shè)恒星成像在水平坐標(biāo)系中的星點(diǎn)坐標(biāo)為（xS，yS，0），則存在方向矢量（xS，yS，-f）滿足下式：

化簡(jiǎn)上式可得：

綜上所述，通過構(gòu)造輔助坐標(biāo)系，依據(jù)傾角儀輸出值和坐標(biāo)變換原理，可得到恒星成像在水平面上的恒星像點(diǎn)坐標(biāo)，并將式（14）稱之為CCD像點(diǎn)修正模型。

4 星點(diǎn)歸算模型

如圖2所示，數(shù)字天頂儀旋轉(zhuǎn)拍攝星圖過程中，數(shù)字天頂儀在不同拍攝位置的視場(chǎng)存在重合，同一顆恒星可能成像在不同拍攝位置的CCD敏感器上。并且在拍攝過程中，CCD平面既有地球自轉(zhuǎn)引起的轉(zhuǎn)動(dòng)，也有數(shù)字天頂儀自身的旋轉(zhuǎn)引起的轉(zhuǎn)動(dòng)，為建立星點(diǎn)歸算模型，需要考慮這兩種轉(zhuǎn)動(dòng)引起的像點(diǎn)CCD坐標(biāo)的變化。

首先考慮地球自轉(zhuǎn)引起CCD平面的轉(zhuǎn)動(dòng)，以水平面為基準(zhǔn)建立CCD北向基準(zhǔn)坐標(biāo)系O-xNyNzN，其中OzN軸過CCD敏感器中心并與垂直軸重合，OxN軸與緯線圈相切，方向指向東，OyN軸與經(jīng)線圈相切，方向指向北。對(duì)于兩個(gè)拍攝位置，假設(shè)數(shù)字天頂儀從前一拍攝位置旋轉(zhuǎn)到下一拍攝位置，CCD北向基準(zhǔn)坐標(biāo)系由 O-xNyNzN變?yōu)?O-x′Ny′Nz′N，如圖2所示。

初值時(shí)刻，水平狀態(tài)下CCD圖像坐標(biāo)系O-xsyszs與CCD北向基準(zhǔn)坐標(biāo)系O-xNyNzN之間存在方位角A，根據(jù)坐標(biāo)轉(zhuǎn)換原理可得：

圖2 CCD北向基準(zhǔn)坐標(biāo)系Fig.2 CCD north reference coordinate

其反變換為：

由第2節(jié)可知，在對(duì)單幅星圖進(jìn)行識(shí)別過程中建立了放射變換模型，其方位角A可由下式得到［8］：

在坐標(biāo)系 O-x′Ny′Nz′N中存在與地球自轉(zhuǎn)軸平行的方向矢量 u=（0，cosδ，sinδ） ，可使坐標(biāo)系 O-xNyNzN繞旋轉(zhuǎn)地軸旋轉(zhuǎn)地球自轉(zhuǎn)角θ后與坐標(biāo)系O-x′Ny′Nz′N平行，則存在四元數(shù) Q：

式中，θ=wret，wre為地球自轉(zhuǎn)角速度。取wre=7.292115×10-5rad/s；t為拍攝任意兩幅星圖間的時(shí)間間隔，可由天頂儀測(cè)得輸出；δ為測(cè)站點(diǎn)的天文緯度。

對(duì)于同一顆恒星，假設(shè)恒星成像在下一位置的CCD北向基準(zhǔn)坐標(biāo)系中的坐標(biāo)為 （x′N，y′N，0） ，由于位于無窮遠(yuǎn)的恒星發(fā)出的星光可認(rèn)為是平行光，則在下一位置的CCD北向坐標(biāo)系O-x′Ny′Nz′N中，像點(diǎn)與焦點(diǎn)構(gòu)成的方向矢量 r=（x′N，y′N， -f）滿足：

由上式可得恒星成像在下一位置CCD北向基準(zhǔn)坐標(biāo)系上的坐標(biāo)為：

坐標(biāo)變換矩陣中包含未知數(shù)δ，由于位置9與位置10拍攝的星圖已經(jīng)單獨(dú)識(shí)別出來，則可得到兩幅星圖上的共有星，故δ角可通過將位置9與位置10之間共有星的CCD圖像坐標(biāo)代入公式（21）求得。

聯(lián)立公式（16）和公式（22），可得：

其次考慮數(shù)字天頂儀鏡筒的旋轉(zhuǎn)，假設(shè)數(shù)字天頂儀旋轉(zhuǎn)角為φ，則：

式（24）稱之為星點(diǎn)歸算模型。

綜上所述，為尋找位置1至位置8拍攝的星圖與位置9拍攝星圖之間的共有星，以及位置11至位置16拍攝的星圖與位置10拍攝的星圖之間的共有星，可按如下步驟進(jìn)行：

（1）將數(shù)字天頂儀拍攝的16幅星圖按式（11）計(jì)算得到輔助坐標(biāo)系下的CCD圖像坐標(biāo)（xF，yF）。

（2）按式算出各位置在水平狀態(tài)下的CCD圖像坐標(biāo) （xS，yS） 。

（3）結(jié)合位置9與位置10得到的星圖上的像點(diǎn)CCD坐標(biāo)按式（24）進(jìn)行計(jì)算，得到轉(zhuǎn)換后的CCD圖像坐標(biāo) （x′S，y′S） 。

（4）若滿足下式，則表示位置9與位置10的星圖與其他位置星圖間存在共有星。

式中，ε為給定的閾值（一般取3pixel左右）。

5 實(shí)驗(yàn)結(jié)果分析

2017年11月，利用數(shù)字天頂儀采用在西安某地拍攝了大量星圖，實(shí)驗(yàn)過程中數(shù)字天頂儀的焦距為（600±4）mm，數(shù)字天頂采用美國(guó)FLI公司生產(chǎn)的ML16803數(shù)碼相機(jī)，CCD敏感器為Kodak公司生產(chǎn)的KAF-16803型面陣CCD。為實(shí)現(xiàn)對(duì)星圖的識(shí)別，首先將依巴谷星表中恒星的位置信息轉(zhuǎn)換拍攝時(shí)刻的視位置，在此基礎(chǔ)上對(duì)星圖進(jìn)行識(shí)別，由于篇幅的限制，這里對(duì)一個(gè)定位循環(huán)的星圖數(shù)據(jù)進(jìn)行分析。其中表1為位置9星圖識(shí)別的結(jié)果，表2為位置1星圖識(shí)別的結(jié)果。

由表1與表2的星圖識(shí)別結(jié)果可知，兩幅星圖之間的共有星為11顆，超過星圖識(shí)別數(shù)量的50%，由此可見星圖間存在著大量的共有星。

結(jié)合公式（24），將位置10星圖上的星點(diǎn)坐標(biāo)歸算到位置16上拍攝的星圖上，為直觀表示出兩幅星圖上的共有星，可得到如圖3所示。

表1 位置9星圖識(shí)別結(jié)果Tab.1 Identification results of star map 9

表2 位置1星圖識(shí)別結(jié)果Tab.2 Identification results of star map 1

圖3 共有星Fig.3 Common stars

由圖3分析可知，位置10與位置16拍攝的星圖上存在14顆共有星，且共有星的分布較為集中。其中將位置10星圖上的共有星按式轉(zhuǎn)換后，與位置16上共有星的CCD圖像坐標(biāo)相比，其中X坐標(biāo)最大相差2.782 pixel，Y坐標(biāo)最大相差2.629 pixel，X坐標(biāo)與Y坐標(biāo)最大差值均不超過3 pixel，由此可知由式（24）構(gòu)造的坐標(biāo)歸算模型滿足尋找星圖間共有星的需求。

為進(jìn)一步說明星圖間存在著共有星，再任意選取一個(gè)定位循環(huán)拍攝的星圖進(jìn)行分析，將位置9星圖上星點(diǎn)坐標(biāo)歸算到順時(shí)針拍攝的8幅星圖上來尋找與星圖1至星圖8與星圖9的共有星，將位置10星圖上的星點(diǎn)歸算到逆時(shí)針拍攝的6幅星圖上來尋找與星圖11至星圖16之間的共有星，如表4所示。

表4 共有星的數(shù)量Tab.4 The number of common star

由表4可知，與星圖8之間的共有星最少15顆，與星圖10之間的共有星最少為14顆。在實(shí)驗(yàn)過程中，利用數(shù)字天頂儀拍攝的單幅星圖識(shí)別的恒星數(shù)量為25顆左右，通過坐標(biāo)歸算模型尋找出星圖間的共有星，再對(duì)剩余的新增星單獨(dú)進(jìn)行識(shí)別，將大大提高星圖識(shí)別的準(zhǔn)確性。例如，假設(shè)一幅星圖上有n顆觀測(cè)星點(diǎn)，則可以組成C3n個(gè)觀測(cè)三角形，有n顆導(dǎo)航星與其匹配，則可以組成C3n個(gè)導(dǎo)航三角形，并且假設(shè)每一個(gè)觀測(cè)三角形需要與每次匹配后剩余的所有導(dǎo)航三角形匹配一次，那么對(duì)星圖上所有的星點(diǎn)組成的觀測(cè)三角形進(jìn)行匹配，需要進(jìn)行次匹配。若采用星點(diǎn)歸算模型找出共有星后，假設(shè)存在j顆共有星，則對(duì)剩余的（n-j）顆新增星進(jìn)行識(shí)別，則需要進(jìn)行次匹配。在本次實(shí)驗(yàn)定位循環(huán)中，位置16提取出21顆觀測(cè)星，與位置10之間存在14顆共有星，若對(duì)所有星點(diǎn)單獨(dú)進(jìn)行識(shí)別需要構(gòu)建1330個(gè)觀測(cè)三角形，并且需進(jìn)行883785次匹配；若提取出14顆共有星后，只需要構(gòu)建35個(gè)觀測(cè)三角形，只進(jìn)行45920次匹配，相比于原來，匹配次數(shù)大大減少，由此可以看出星圖識(shí)別的速度得到較大的提高。

6 結(jié) 論

通過建立星點(diǎn)歸算模型，能夠依據(jù)位置9與位置10拍攝星圖上星點(diǎn)信息找出其他星圖上的共有星，單獨(dú)對(duì)其他星圖的新增星進(jìn)行識(shí)別，將大大減少匹配三角形的個(gè)數(shù)和匹配的次數(shù)，提高星圖識(shí)別的速度。