崔平遠(yuǎn)，賈 賀，朱圣英

(1. 北京理工大學(xué)深空探測技術(shù)研究所，北京 100081；2. 深空自主導(dǎo)航與控制工信部重點(diǎn)實(shí)驗(yàn)室，北京 100081；3. 飛行器動(dòng)力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京 100081)

0 引 言

小天體探測具有重要的科學(xué)價(jià)值，是近些年世界各國深空探測的熱點(diǎn)問題。自主導(dǎo)航與控制系統(tǒng)是成功實(shí)現(xiàn)小天體近距離探測任務(wù)的前提和保證[1-2]。光學(xué)導(dǎo)航(OPNAV)憑借自主性強(qiáng)、精度高等優(yōu)點(diǎn)，在深空探測自主導(dǎo)航方面已經(jīng)有了廣泛應(yīng)用[3-4]。小天體探測一般通過拍攝目標(biāo)天體圖像，然后提取目標(biāo)天體的形心或光心，以形心和光心代替質(zhì)心得到相機(jī)至天體質(zhì)心的視線(LOS)信息進(jìn)行導(dǎo)航[5-6]。不規(guī)則小天體不同的成像特點(diǎn)對小天體形心/光心提取精度有很大影響，而小天體形心/光心提取的精度直接影響到自主光學(xué)導(dǎo)航的精度，所以高精度的形心/光心提取算法是至關(guān)重要的。

針對形心提取問題，Christian等[7-8]提出了一種目標(biāo)天體形心和視直徑(CAD)的圖像處理算法，首先利用橢圓和圓擬合算法求取目標(biāo)形心，得到視線方向，然后利用目標(biāo)天體的形狀模型和投影圖像進(jìn)行比對得到視直徑信息。但是提前獲得小天體形狀和尺寸的精確模型是很困難的，另外由于小天體形狀不規(guī)則，導(dǎo)致在相機(jī)像平面的投影不規(guī)則，所以很難將投影圖和實(shí)際模型進(jìn)行對比分析?；谔祗w邊緣進(jìn)行橢圓擬合是實(shí)現(xiàn)形心提取的主要方法。常用的橢圓擬合方法主要有3類，一類是基于HOUGH變換的橢圓擬合方法[9]，另一類方法是基于不變矩的方法[10]，再一類則是基于最小二乘的擬合方法[11]。橢圓擬合方法僅適用于真實(shí)天體邊緣近似規(guī)則橢圓的情況，當(dāng)天體邊緣不規(guī)則時(shí)，則會(huì)存在較大誤差，造成橢圓擬合的誤匹配，從而得到誤差較大的天體形心位置坐標(biāo)。

針對光心提取問題。高斯擬合法及其改進(jìn)算法是近些年較為常用的光心提取方法[12-13]，包括基于不飽和點(diǎn)的高斯曲面擬合法[14]，利用一階導(dǎo)數(shù)零交叉點(diǎn)的方法確定峰值像素坐標(biāo)分布后利用高斯擬合法求解等[15]。在理想情況下，光斑的強(qiáng)度分布應(yīng)該符合高斯分布，所以高斯曲面擬合法提取精確度高。然而當(dāng)亮點(diǎn)的強(qiáng)度分布為非高斯的情況則精度較差，此時(shí)只能選用灰度重心法提取小天體的光心。該方法通過在一個(gè)預(yù)先確定的圖像子區(qū)域內(nèi)提取超過亮度閾值的所有像素，來確定該區(qū)域的亮度中心，該方法簡單快速，但提取精度有限。

針對以上問題，考慮到不規(guī)則小天體不同的成像特點(diǎn)對圖像處理算法的提取精度影響很大，因此需要對小天體的不規(guī)則程度進(jìn)行判斷，評價(jià)擬合精度，分布特性等指標(biāo)，進(jìn)而選取最合適的擬合方法。本文提出了一種小天體不規(guī)則度定義方法。小天體不規(guī)則度可以表示小天體的不規(guī)則程度，同時(shí)可以作為小天體形心/光心提取方法的選取指標(biāo)，實(shí)現(xiàn)對不同成像特點(diǎn)的算法確定。另外，針對小天體形狀不規(guī)則度超過設(shè)定閾值的情況(即小天體為不規(guī)則小天體)，本文提出了一種基于最小外接圖形的形心提取方法，實(shí)現(xiàn)不規(guī)則小天體形心精確提取。該方法受小天體自旋影響小，能夠提高小天體形心的提取精度，進(jìn)而提高自主光學(xué)導(dǎo)航的精度。

本文的整體布局如下：首先提出了小天體不規(guī)則度的概念，并給出了定義。接著針對小天體探測小天體成像特點(diǎn)，給出了基于小天體不規(guī)則度的算法選取準(zhǔn)則。針對不規(guī)則度較大的小天體，提出了一種基于最小外接圖形的小天體形心提取方法，并對該方法進(jìn)行了詳細(xì)說明；然后對本文提出的圖像處理方法及選取準(zhǔn)則進(jìn)行數(shù)值仿真驗(yàn)證和分析；最后給出了文章結(jié)論。

1 小天體不規(guī)則度定義

不規(guī)則圖形即沒有規(guī)則外形的圖形，其面積周長等參數(shù)的計(jì)算沒有規(guī)律可循。亮度的不規(guī)則體現(xiàn)在亮斑每一個(gè)像素點(diǎn)灰度值分布的非高斯性。針對距離小天體遠(yuǎn)近兩種不同的成像特點(diǎn)，將小天體形狀不規(guī)則度和亮度不規(guī)則度相結(jié)合，提出了一種小天體不規(guī)則度(Asteroid irregularity degree，AID)的定義，并作為小天體質(zhì)心提取的選取指標(biāo)，小天體不規(guī)則度的表達(dá)式如下所示：

(1)

1.1 小天體形狀不規(guī)則度

對于不規(guī)則的小天體，因?yàn)檫吘墳椴灰?guī)則曲線，因此無法直接計(jì)算邊緣點(diǎn)的曲率及投影圖形的每一個(gè)內(nèi)角。因此，本文提出了一種基于凸包的小天體形狀不規(guī)則度的定義方法。根據(jù)小天體形狀偏離對稱性的程度表征形狀不規(guī)則度。首先對相機(jī)拍攝到的小天體原始圖像進(jìn)行閾值分割，接著計(jì)算目標(biāo)小天體的凸包，然后根據(jù)小天體的凸包，利用綜合判別因素來描述小天體形狀的不規(guī)則程度。其中小天體形狀不規(guī)則度的定義如下：

(2)

以三角形為例，首先給出軸對稱判別系數(shù)Rs的定義。如圖1所示，黑色粗箭頭直線代表經(jīng)過三角形ABC重心O點(diǎn)的對稱軸MN，三角形ABC為平面投影圖形，三角形A1B1C1為ABC的軸對稱圖形，S為投影三角形ABC的面積，s是投影與其軸對稱圖形重合部分(陰影區(qū))的面積。通過旋轉(zhuǎn)對稱軸MN可以得到陰影區(qū)面積的最大值smax，進(jìn)而可以得到投影圖形的軸對稱判別系數(shù)Rs，表達(dá)式如下：

圖1 軸對稱性判別系數(shù)Fig.1 Schematic diagram of ellipse parameters

Rs=1-smax/Sa

(3)

對于尺寸不規(guī)則判別系數(shù)Rd，其表達(dá)式如下：

(4)

式中:di為凸包第i個(gè)頂點(diǎn)到等效中心的距離，等效中心可以通過計(jì)算凸包頂點(diǎn)的一階矩得到。d0是由當(dāng)前圖像所求得凸包的等效直徑。

對于形狀不規(guī)則判別系數(shù)Ra，其表達(dá)式如下：

(5)

式中:θi是凸包多邊形的第i個(gè)內(nèi)角，θ0為凸包對應(yīng)的相同邊數(shù)的正多邊形的內(nèi)角。

對小天體邊緣信息進(jìn)行提取后，在所建立的初步匹配特征中可能存在誤匹配[16]。誤匹配邊緣曲線的存在極大地影響探測器運(yùn)動(dòng)參數(shù)的求取，使估計(jì)結(jié)果出現(xiàn)較大偏差。下面以橢圓為例對擬合殘差系數(shù)進(jìn)行說明。這里給出橢圓擬合殘差定義：

(6)

當(dāng)橢圓擬合的殘差足夠小時(shí)(低于閾值下限)，邊緣點(diǎn)和橢圓輪廓的近似程度高，采用基于Hough變換的橢圓擬合方法具有很高的可靠性。當(dāng)橢圓擬合殘差大于閾值下限時(shí)，邊緣點(diǎn)與橢圓的輪廓差距較大，即采用橢圓擬合方法的可靠性較差，此時(shí)需要采用適用于不規(guī)則邊緣的處理算法進(jìn)行目標(biāo)天體形心提取。同樣可以給出矩形模板和三角形模板的擬合殘差。

(7)

因此，針對小天體二維圖像，小天體形狀不規(guī)則度的具體表達(dá)形式如下：

(8)

針對三維的小天體，小天體不規(guī)則度的具體表達(dá)形式如下：

(9)

1.2 小天體亮度不規(guī)則度

距離小天體較遠(yuǎn)時(shí)小天體在相機(jī)中的成像為一個(gè)亮點(diǎn)或亮斑，一般其灰度值成高斯分布。本節(jié)提出了一種亮度不規(guī)則度來描述小天體成像亮度分布的不規(guī)則程度。其中小天體亮度不規(guī)則度的定義如下：

(10)

偏斜度是對統(tǒng)計(jì)數(shù)據(jù)分布偏斜方向及程度的度量[17]。頻數(shù)分布分為對稱和非對稱(偏態(tài))兩種。頻數(shù)分布的示意圖如圖2所示。在偏態(tài)分布中，當(dāng)偏斜度為正值時(shí)，分布不對稱，眾數(shù)位于算術(shù)平均數(shù)的左側(cè)或右側(cè)。峭度是一個(gè)無量綱參數(shù)或統(tǒng)計(jì)量，用于量化信號相對于高斯分布的分布形狀[18]。信號分布的形狀分為尖銳，平坦和高斯分布。當(dāng)信號形狀與高斯分布相同時(shí)，峭度為零；當(dāng)信號形狀比高斯分布更尖銳時(shí)，峭度為正；當(dāng)信號形狀比高斯分布更平坦時(shí)，峭度為負(fù)。X方向上的偏斜度系數(shù)和峭度系數(shù)表達(dá)式如下。

圖2 圖像灰度頻率分布圖Fig.2 Geometric sketch of frequency distribution of image gray values

(11)

因此，小天體亮度不規(guī)則度表達(dá)式如下：

(12)

式中:Rcy，Rky分別為Y方向上偏斜度系數(shù)和峭度系數(shù)。根據(jù)1.1節(jié)和1.2節(jié)給出的小天體形狀不規(guī)則度和小天體亮度不規(guī)則度的定義，可以計(jì)算出小天體不規(guī)則度，具體表達(dá)式如下：

μ2((Rcx+Rcy)+(Rkx+Rky))

(13)

當(dāng)目標(biāo)小天體圖像為一個(gè)亮點(diǎn)或亮斑時(shí)，μ1=0；反之μ1=1。根據(jù)具體任務(wù)設(shè)計(jì)與目標(biāo)天體參數(shù)可以判斷目標(biāo)天體成像特性。本文中假定當(dāng)小天體所占像素個(gè)數(shù)小于50時(shí)即為亮斑。小天體不規(guī)則度可以評定小天體的不規(guī)則程度，并可以作為圖像處理算法的選取指標(biāo)。

2 基于小天體不規(guī)則度的質(zhì)心提取算法

小天體探測一般很難提前獲取目標(biāo)小天體的精確形狀模型?？紤]小天體不規(guī)則度的適用性及形心提取方法的可行性，以小天體接近段為例，本節(jié)提出了一種基于小天體不規(guī)則度的算法選取準(zhǔn)則，另外針對小天體不規(guī)則度較大的一類小天體，提出了一種最小外接圖形的質(zhì)心提取方法。

2.1 基于小天體不規(guī)則度的算法選取

當(dāng)探測器距離目標(biāo)小天體相對較近時(shí)，小天體在像平面中的投影為一個(gè)二維圖形，當(dāng)真實(shí)的天體邊緣為規(guī)則邊緣時(shí)，采用基于Hough檢測提取形心的算法具有較高的精度，而當(dāng)天體邊緣不規(guī)則時(shí)則具有較大提取誤差。當(dāng)探測器距離目標(biāo)小天體較遠(yuǎn)時(shí)，小天體在像平面中的成像為一個(gè)亮點(diǎn)或亮斑。理想情況下，亮點(diǎn)的灰度值分布應(yīng)該符合高斯分布，所以采用高斯曲面擬合法有較高的精度，而當(dāng)亮點(diǎn)灰度值不符合高斯分布時(shí)則具有較大提取誤差。因此，本文根據(jù)相機(jī)拍攝的不同時(shí)刻的不規(guī)則小天體圖像，提出了一種基于小天體不規(guī)則度的圖像處理算法選取準(zhǔn)則，選取最合適的小天體形心/光心提取算法。算法選取流程如圖3所示。

圖3 選取準(zhǔn)則Fig.3 Selection criteria based on asteroid irregularity degree

閾值是區(qū)分目標(biāo)規(guī)則和不規(guī)則的臨界值。它是通過基于多場景和多任務(wù)圖像處理結(jié)果的經(jīng)驗(yàn)獲得的。在本文中，根據(jù)作者的經(jīng)驗(yàn)，將火星遠(yuǎn)，中，近序列圖像計(jì)算出的不規(guī)則度平均值的兩倍用作閾值，即4.27。應(yīng)該指出的是，該值不一定是最優(yōu)值，閾值選擇仍需要進(jìn)一步的研究。

2.2 不規(guī)則小天體形心提取算法

針對小天體不規(guī)則度較大的一類小天體，本文提出了一種基于最小外接圖形的形心提取方法，該方法先求投影圖像的邊緣，再根據(jù)邊緣點(diǎn)求出目標(biāo)小天體的凸包頂點(diǎn)，然后根據(jù)小天體形狀不規(guī)則度的數(shù)值及擬合殘差計(jì)算凸包頂點(diǎn)的最小外接圖形提取形心。該方法確保了計(jì)算結(jié)果的精確性，提高了不規(guī)則小天體形心的提取精度，下面以圖4(a)為原始圖像，對該提取方法進(jìn)行詳細(xì)說明：

1) 圖像預(yù)處理

圖像預(yù)處理技術(shù)的目的是使相機(jī)拍攝的圖像更適合于后續(xù)處理，一般通過增強(qiáng)感興趣的區(qū)域同時(shí)抑制沒有實(shí)際作用的區(qū)域?qū)崿F(xiàn)。本文采用中值濾波器進(jìn)行導(dǎo)航圖像的預(yù)處理。

2) 邊緣檢測

邊緣檢測就是利用離散化函數(shù)根據(jù)二維灰度矩陣梯度向量來尋找圖像灰度矩陣的灰度躍變位置，然后在圖像中將這些點(diǎn)連接起來就構(gòu)成了目標(biāo)的圖像邊緣。Canny邊緣檢測是一種常見的邊緣提取方法[19]，它利用特定的數(shù)字掩碼對圖像進(jìn)行濾波運(yùn)算，具有提取速度快的特點(diǎn)，本文利用Canny算子提取目標(biāo)的邊緣。邊緣檢測的結(jié)果如下圖4(b)所示。

3) 計(jì)算凸包

凸包是指覆蓋平面坐標(biāo)系內(nèi)若干點(diǎn)的面積最小的凸多邊形[20]，根據(jù)步驟2提取的圖像中目標(biāo)小天體邊緣輪廓，采用格雷厄姆掃描Graham-Scan算法進(jìn)行小天體圖像邊緣的凸包計(jì)算，該算法的時(shí)間復(fù)雜度為O(n)，小天體凸包計(jì)算結(jié)果如圖4(c)所示。

圖4 處理結(jié)果Fig.4 Processing results

4) 計(jì)算最小外接圖形提取形心

根據(jù)計(jì)算得到的小天體凸包，根據(jù)擬合殘差計(jì)算凸包的最小外接圖形，最小外接圖形的形心即為目標(biāo)小天體的質(zhì)心，最終實(shí)現(xiàn)導(dǎo)航信息的提取。根據(jù)得到的小天體凸包，利用哈奇揚(yáng)Khachiyan算法計(jì)算凸包頂點(diǎn)的最小外接橢圓(MEE)[21]。通過計(jì)算面積最小的外接矩形即可得到最小外接矩形(MER)。通過O′Rourke’s算法可以求解凸多邊形的最小外接三角形(MET)[22]。該算法通過for循環(huán)使某條邊依次與多邊形的每個(gè)邊齊平，而不管該邊是否是下一個(gè)要接觸的邊。因此，該算法搜索局部極小點(diǎn)的超集，即平面三角形。

3 仿真分析

本節(jié)將使用深空導(dǎo)航光學(xué)成像半物理仿真系統(tǒng)對本文提出的基于小天體不規(guī)則度的質(zhì)心提取方法進(jìn)行仿真驗(yàn)證。相機(jī)的分辨率為1024×1024。為了解決目標(biāo)天體建模問題，選用3 ds max進(jìn)行小天體建模工作。通過對3 ds格式的三維模型進(jìn)行修改與貼圖輸出lwo格式，再通過LwConvert格式轉(zhuǎn)換工具將模型格式進(jìn)行轉(zhuǎn)換，最終通過Open GL調(diào)用三維數(shù)字模型得到仿真模擬圖像。其中小天體模型如圖5所示。不失一般性，將模擬小天體質(zhì)心的真實(shí)位置設(shè)置在圖像的中心，以便更直觀地比較處理結(jié)果。通過將真實(shí)位置與所提算法的提取結(jié)果進(jìn)行比較，可以獲得提取誤差。另外，對于不規(guī)則的小天體圖像，對基于小天體不規(guī)則度的質(zhì)心提取方法進(jìn)行了仿真分析。

圖5 三維仿真模型Fig.5 Three-dimensional simulation model

3.1 小天體成像為一個(gè)亮點(diǎn)

圖6 光點(diǎn)仿真圖像Fig.6 Synthetic simulated images of light spot

表1 光心提取結(jié)果對比Table 1 Comparison of optical center extraction results

通過仿真可以看出，小天體不規(guī)則度能夠評價(jià)小天體光點(diǎn)亮度分布的不規(guī)則程度，并作為算法選取的評價(jià)指標(biāo)。對于不同的成像情況，圖像處理精度均小于0.35個(gè)像素，證明本文提出的基于小天體不規(guī)則度的光心提取方法能夠?qū)崿F(xiàn)不規(guī)則小天體光心的精確提取。

3.2 小天體成像為不規(guī)則圖形

為了驗(yàn)證選擇標(biāo)準(zhǔn)的可靠性和所提算法的提取精度，對于圖4 (a)中的小天體圖像，圖7中提供了不同光照條件和小天體自旋狀態(tài)的仿真圖像。根據(jù)選擇標(biāo)準(zhǔn)選擇適當(dāng)?shù)膱D像處理算法，結(jié)果如圖8和表2所示。圖8(f)橢圓表示EHT算法的擬合結(jié)果，而圖8(a～e)橢圓和矩形表示基于MEE和MER算法的處理結(jié)果。十字表示小天體的真實(shí)位置，星號表示圖像處理的質(zhì)心提取結(jié)果。

通過圖8和表2可以看出基于小天體不規(guī)則度選取后的提取算法有很高的提取精度，形心提取誤差小于1個(gè)像素，說明了選取準(zhǔn)則的有效性，同時(shí)證明本文提出的基于外接圖形的形心提取算法對于不規(guī)則度較大的一類小天體有較高的提取精度。

為了進(jìn)一步說明選擇方法的有效性，圖7中的小天體圖像使用了四種不同的圖像處理算法。圖9中顯示了圖像處理誤差曲線。

圖7 小天體不同狀態(tài)模擬圖像Fig.7 Synthetic images of asteroid with different states

圖8和圖9的仿真結(jié)果表明，對于不同的光照和旋轉(zhuǎn)狀態(tài)，所選擇的算法具有最高的提取精度，這證明小天體不規(guī)則度可以有效地選擇不同條件下的圖像處理算法。根據(jù)表2的處理結(jié)果，圖像處理精度基本上優(yōu)于1個(gè)像素，并且提取誤差的最大值為1.0537。為了進(jìn)一步驗(yàn)證提取算法在自主光學(xué)導(dǎo)航中的可行性，將圖像處理誤差的最大值作為相機(jī)測量誤差，并基于圖像處理結(jié)果對導(dǎo)航性能進(jìn)行了驗(yàn)證和分析。

圖8 圖像處理結(jié)果Fig.8 Image processing results

圖9 圖像處理結(jié)果對比Fig.9 Comparison of different image processing algorithm

表2 形心提取結(jié)果對比Table 2 Comparison of shape center extraction results

假設(shè)攝像機(jī)分辨率為1024×1024，視場為1.5°，小天體固定坐標(biāo)系中探測器的初始狀態(tài)為[45,000 km，45,000 km，45,000 km，-2 km/s，-1 km/s，-1 km/s]，探測器的初始位置誤差為[30，30，30] km，探測器的初始速度誤差為[1，1，1] m/s，濾波周期為30 s，仿真時(shí)間為1.5×104s。根據(jù)小天體接近段的動(dòng)力學(xué)模型和相機(jī)測量模型，探測器的估計(jì)誤差如圖10所示。

圖10 導(dǎo)航狀態(tài)估計(jì)誤差曲線Fig.10 State estimation error curve

從仿真結(jié)果分析可知:

1) 針對不規(guī)則度較大的小天體，本文提出的基于邊緣凸包的最小外接圖形進(jìn)行形心提取的方法具有較高的提取精度，提取誤差小于1個(gè)像素，同時(shí)受光照影響小。

2) 基于Hough變換的形心提取方法對于規(guī)則圖形邊緣具有較高的提取精度，而對于不規(guī)則圖形邊緣精度較差。本文提出的基于小天體不規(guī)則度的圖像處理方法選取準(zhǔn)則，能夠選擇最優(yōu)的形心提取方法。通過仿真可以看出，小天體不規(guī)則度可以評價(jià)小天體的不規(guī)則程度，實(shí)現(xiàn)算法的選取。對于不同的成像情況，選取后的圖像處理精度最高，證明本文提出的基于小天體不規(guī)則度可以實(shí)現(xiàn)不同算法的有效選取。

3)導(dǎo)航誤差曲線表明導(dǎo)航系統(tǒng)具有良好的收斂性。接近段探測器位置估計(jì)誤差小于1.19 km，速度估計(jì)誤差小于0.12 m/s。因此，本文提出的質(zhì)心提取算法的精度滿足小天體接近段導(dǎo)航信息提取的精度要求。

4 結(jié) 論

本文針對小天體形心、光心提取的問題，提出了小天體不規(guī)則度的定義，能夠表征小天體的不規(guī)則程度。小天體不規(guī)則度可以作為圖像處理算法的選取指標(biāo)，針對不同的成像特點(diǎn)采用不同的處理方法，實(shí)現(xiàn)對導(dǎo)航形心或光心位置的精確提取。針對小天體不規(guī)則度較大的一類小天體，本文提出了一種基于最小外接圖形的形心提取方法，該方法受光照和天體自旋影響小，提高了不規(guī)則小天體形心的提取精度，最終提高了導(dǎo)航精度。