先驗信息輔助的運動模糊星圖星點提取方法

何貽洋，王宏力，馮 磊，由四海，劉 珂，高迎東

(火箭軍工程大學 導彈工程學院，西安 710025)

星敏感器是目前公認的精度最高的姿態敏感器之一[1].在航空、航天領域，無論是近地衛星或者深空探測，還是彈道導彈或者小衛星上，星敏感器都得到了廣泛的應用.星敏感器通常捷聯安裝于上述的載體上，若載體受到外界環境的沖擊或者進行大姿態角機動時，星敏感器拍攝的星圖會出現運動模糊現象，導致星圖信噪比降低，可能嚴重影響星點提取成功率和精度，進而降低星圖識別成功率[2-3].

目前對運動模糊星圖處理的主流思路是先復原星圖再進行后續的星點提取、識別等過程.根據上述思路，運動模糊星圖處理的核心在于星圖復原，故而有不少學者開展了星圖復原的研究.目前星圖復原方法主要是非盲圖像復原方法，非盲圖像復原方法應用最為廣泛的是逆濾波方法[4-5]、維納濾波方法[6-7]和Lucy-Richardson(L-R)方法[8-9].逆濾波方法是應用最早的圖像復原方法，但該方法缺點很明顯，即對圖像噪聲非常敏感；維納濾波方法可以實現最小均方復原效果，但該方法需已知圖像噪聲特性的峰值信噪比較高時才能獲得較好的復原效果；L-R方法是一種統計學方法，如果圖像的先驗概率模型和條件概率模型與實際不符時，圖像復原效果會受到極大影響.另外，非盲圖像復原方法受圖像先驗信息的準確性影響很大，故而對圖像的模糊核函數估計也是研究的重點，Radon變換方法[10-11]、基于星像運動軌跡的估計方法[12-14]是當前模糊核函數估計的常用方法.傳統的Radon變換方法的估計精度受星圖噪聲影響很大，Jiang等[15]提出了改進的Radon變換方法，該方法基于模糊星圖噪聲特性并且結合了Z函數和雙線性插值方法，取得了較好的估計效果；吳小娟等[12]基于繞星敏感器光軸或者非光軸轉動時的星像運動軌跡模型進而得到相應條件下星圖模糊核函數的解析表達式；廖育富等[14]根據星像運動軌跡模型提出了一種操作簡單利于工程化實現的模糊核函數估計方法，但估計精度不高且容易受到噪聲影響.上述都是基于星圖自身信息估計模糊核函數，但星圖本身的先驗信息有限，所以近幾年有學者提出基于慣性導航系統(inertial navigation system，INS)輔助對模糊核函數進行估計[16-18]，取得了較好的估計效果.

綜上所述，對于運動模糊星圖的研究聚焦于模糊核函數的估計和復原方法，期望獲得質量較好的復原星圖以保證星點提取精度和星圖識別成功率，但星圖復原過程會占據星敏感器一定的計算資源，影響星敏感器的數據更新率.因此，本文提出了一種基于INS和星敏感器的先驗信息直接提取運動模糊星圖上星點的方法.該方法首先利用INS和星敏感器的姿態先驗信息確定局部區域生長的“種子”點，從而得到運動模糊星圖上的星像區域；然后基于先驗信息建立了星點提取誤差修正模型，補償利用質心法直接提取運動模糊星圖星點的坐標誤差.該方法可以跳過星圖復原操作直接對運動模糊星圖進行星點提取，這對于保證高動態環境下星敏感器的數據更新率具有重要意義.

1 方法分析

高動態條件下星敏感器成像受到影響，星圖會出現模糊現象，圖1是模擬的一幅運動模糊星圖，模擬星圖的灰階為8位，灰度范圍為0～255，添加的高斯白噪聲均值為0，均方差為6(灰度值)，載體的角運動速度設置為ωx=-2°/s，ωy=-3°/s，ωz=2°/s.分別利用L-R方法[8-9]和維納濾波方法[6-7]進行復原處理后，得到局部放大的復原星圖以及相應的灰度分布如圖2所示.

從上述星圖復原的結果來看，利用維納濾波方法和L-R方法復原星圖后，并不能將星點恢復到靜態條件下的高斯彌散斑，但星點“拖尾”軌跡相對縮短，利用質心法對復原星圖上星點提取的誤差小于直接對原始星圖上星點提取的誤差.如果能夠估計出直接利用質心法提取模糊星圖上星點的誤差，然后對星點坐標進行補償修正，這就能跳過星圖復原這個過程，從而大大節省了時間，保證星敏感器數據更新率，同時仍能獲得較好的星點提取精度.基于上述思想，本文提出了基于INS輔助的直接對運動模糊星圖星點提取的方法.

圖2 不同方法的星圖復原結果局部放大和灰度分布圖

星敏感器常與INS組合使用，雖然INS存在姿態輸出誤差隨時間累積的缺點，但其短時間的姿態輸出精度很高，文獻[16-18]中將INS的姿態輸出引入到模糊核函數估計中取得了較好的效果，本文進一步利用INS提供的姿態先驗信息，提出一種直接對運動模糊星圖進行星點提取的方法，具體思路如下：

1)基于INS的短時姿態輸出和星敏感器在上一幀星圖的姿態先驗信息，預測出恒星在當前幀星圖的粗略星點坐標，將其作為種子點進行局部區域生長，從而獲得星像區域；

2)利用質心法[3]直接在星像區域計算星點質心坐標；

3)根據INS得到在曝光時間內星敏感器的轉動角速度，計算星圖上星像的運動軌跡長度，根據星像運動軌跡與直接質心法的星點提取誤差之間的關系，修正星點坐標. 本文方法的實施步驟如圖3所示.

1.1 基于先驗信息的局部區域生長

局部區域生長法[19]的三要素是種子點、生長準則和終止條件，而種子點是其中最為關鍵的一環.利用局部區域生長可以快速而準確獲取運動模糊星圖上各個星像的范圍，但必須首先解決如何獲取滿足條件的種子點的問題.單幀星圖的先驗信息有限，本文基于INS從t0→t0+Δt時刻的姿態輸出并結合t0時刻星敏感器輸出的姿態信息，得到t0+Δt時刻的粗略姿態，解得t0+Δt時刻星敏感器光軸指向，進而預測出t0+Δt時刻恒星在成像平面上的粗略坐標，將其作為種子點，利用局部區域生長法獲得運動模糊星圖上星像范圍.為質心法提取星點坐標奠定基礎，具體思路如下：

假設星敏感器根據當前這一幀星圖獲得t0時刻的姿態四元數為q，利用INS的姿態輸出得到從t0→t0+Δt的姿態機動四元數為q′，那么根據四元數的乘法順序和規則，t0+Δt時刻星敏感器的姿態四元數可以由下式得到:

q″=q′q.

(1)

文獻[20]中關于姿態四元數與星敏感器光軸指向(A,D)和星敏感器的滾動角θ之間的關系式：

(2)

(3)

式中:A∈[0,2π)，D∈[-π/2,π/2]，θ∈[0,2π).

得到星敏感器光軸指向(A,D)之后，預測星敏感器在t0+Δt時刻星圖上星像坐標就是建立恒星在星圖上質心坐標與在天球坐標系中的坐標(α,δ)之間的數學模型，利用t0時刻的星敏感器的先驗信息可以得到(α,δ).根據星敏感器成像原理以及星敏感器坐標系與天球坐標系之間的坐標變換關系，在星敏感器光軸指向(A,D)、恒星在天球坐標系中的坐標(α,δ)和星敏感器視場確定的情況下，理論上出現在星圖上的恒星位置也隨之確定.得到此時恒星的星像坐標為

(4)

如果考慮到在實際中星敏感器的滾動角θ一般不為0，那么則有

(5)

根據式(1)～(5)，可以預測得到t0+Δt時刻星圖上星像坐標，將這些星像坐標記為{Si}，其中Si=(xi,yi)，i=1,…,N，N為星圖上星像個數. 以{Si}作為局部區域生長的種子點，從而獲得模糊星圖上各個星像的實際范圍關于區域生長法的原理可具體詳見文獻[19]，在此不再贅述.

1.2 星點提取誤差修正模型

質心法具有操作簡單、魯棒性好的優點，是目前常用的星點提取方法.星像的能量分布對質心法的精度有很大影響，根據文獻[3]，星像能量分布服從二維高斯分布時，質心法精度較高.在高動態環境下，星像不再是靜態時的高斯彌散斑，而是形成了一條“拖尾”.在曝光時間內Te內，星像區域的灰度也不再服從二維高斯分布，因此直接用質心法提取星點的誤差會大大增加.為了對該誤差進行修正，根據文獻[21]的研究，首先應對高動態條件下的運動模糊星圖的灰度分布進行分析建模，在此基礎上分析推導直接用質心法提取星點的誤差.

將星敏感器的曝光時間Te等分為N個微小的區間[ti,ti+1]，i=0,1,…,N-1，Δt=ti+1-ti=Te/N，t0=0，tN=Te，由于[ti,ti+1]非常小，因而可近似認為在[ti,ti+1]內，恒星在成像平面上所形成星像的灰度仍服從二維高斯分布，設在[ti,ti+1]的二維高斯分布的中心點坐標為(xi,yi)，在[ti,ti+1]內，恒星在成像平面上所形成星像的灰度分布如下：

式中:Gi=G0/N，G0是與曝光時間Te內星像總能量相關的常值；Gi(x,y)為[ti,ti+1]內成像平面上的坐標(x,y)處的灰度值;σ為高斯彌散斑半徑，代表了彌散程度，根據二維高斯分布的特點以及3σ準則，星像99%以上的能量集中于距離星點的3σ范圍以內.

假設各個高斯分布的彌散斑半徑相等，將各個微小曝光時間內形成的星像灰度分布函數進行疊加，得到在曝光時間Te內恒星在成像平面上星像的灰度分布函數如下：

利用質心法可以得到[ti,ti+1]內對應的中心點(xi,yi)的計算公式如下：

將式(6)、(7)分別代入式(8)、(9)，則有

根據文獻[14,21-22]研究可知繞光軸的角速度對星圖的影響很小，因而星圖上星像運動軌跡近似為勻速直線運動上的一段，假設星像在成像平面的OX及OY軸上的運動速度分別為vx≈-fωy/dCCD，vy≈fωx/dCCD[12]，ωx、ωy為星敏感器轉動角速度大小，逆時針轉動為正.則在曝光時間Te內，星像運動軌跡在在成像平面的OX及OY軸上的分量為Lx=vxTe，Ly=vyTe，且有

式中i=0,1,…,N-1，Δt=ti+1-ti.

同理，可得：

(10)

可見，星圖運動模糊引起的星點坐標誤差等于星像運動軌跡長度的一半.利用INS可得到星敏感器轉動的角速度，進而估算出星像運動軌跡長度，根據式(10)可對星點坐標進行誤差補償，從而獲得較為準確的星點坐標.

2 方法評價

對于運動模糊星圖處理問題，本文提出一種基于INS輔助直接對運動模糊星圖進行星點提取的方法.INS的輔助作用主要體現在以下兩個方面：1)根據INS的姿態輸出可以得到相鄰兩幀星圖之間星敏感器的姿態角變化量，結合星敏感器在上一幀星圖的姿態輸出，可以預測出當前這幀星圖上星點的粗略坐標；2)基于INS的角速度輸出可估算曝光時間內的星像運動軌跡長度.

本文方法需要利用恒星在上一幀星圖的質心坐標和INS系統的角速度輸出等先驗信息來預測恒星在當前幀星圖中的坐標，作為局部區域生長所需的“種子”點.在實際環境中，先驗信息也不可避免的存在誤差，那么預測的星點坐標也會產生誤差，但若“種子”點仍在實際星像區域范圍內，那么局部區域生長就能準確的獲得實際星像范圍.由此可見，局部區域生長法對“種子”點的精度要求不高，這在一定程度上增強了本文方法對先驗信息誤差的容錯能力；局部區域生長只需關心存在星像的局部星圖，星圖其他區域不需做進一步處理，這相當于對星圖做了一次去噪處理；本文方法避免了掃描整幅星圖，而且不需對運動模糊星圖進行復原操作，相比于先復原星圖再提取星點的方法，本文方法能夠大大節省星點提取時間，故本文方法還具有快速性好、抗噪能力強的優點.

值得注意的是本文方法要求已知恒星的先驗信息，也就是需要已知恒星在上一幀星圖的星點坐標，對應的星敏感器工作模式是星跟蹤模式.如果星敏感器處于星捕獲模式下，也就是有新的恒星進入視場，而新進入視場內的恒星難以獲得其在上一幀星圖中的先驗信息，那么這種情況下本文方法可能將不再適用.

3 仿真實驗及分析

3.1 仿真條件設置

采用變形星敏感器視場設置為9°×9°，焦距為65.76 mm，像元大小為20 um，面陣大小為512×512，曝光時間為100 ms；從SKY2000主星表篩選出星等于小于6的恒星，并剔除其中的雙星與變星，剩余4 908顆組成本文的星表；星敏感器角速度設置為(-2°/s -3°/s 2°/s)，角速度方向記順時針轉動為正；模擬星圖的灰階為8位，灰度范圍為0～255，添加的高斯白噪聲的均值為0，均方差為6(灰度值)；仿真實驗工具是MATLAB(R2015a版)，仿真實驗的電腦配置為2.5 GHz Core(TM)i7-6500U，8.0 GB RAM.

利用Monte-Carlo方法隨機產生100組星敏感器在t0時刻的姿態四元數{qi}和100組t1時刻的光軸指向{(αi,δi)}i=1,…,100，從t0→t1時刻的星敏感器姿態變化四元數{q′i}可利用INS系統的輸出而得到，根據式(2)～(5)可得出t1時刻星敏感器成像平面上星像的粗略坐標，將此作為本文方法中局部區域生長法的種子點.

3.2 仿真實驗分析

3.2.1 仿真實驗1

為了考察“種子”點的誤差對于本文方法結果的影響大小，也即考察本文方法對先驗信息誤差的容錯能力，設計了本仿真實驗.

為了便于對一幅星圖上的星點提取精度進行描述，以星圖上提取到的星點坐標與相應的理想星點坐標距離的平均值作為衡量星點提取精度高低的指標，定義如下：

可以假設“種子”點的誤差是服從均值為0的高斯白噪聲，共設置10組不同大小的“種子”點的誤差，分別為0.5,1.0,…,5.0像素.根據t1時刻的光軸指向{(αi,δi)}得到100幅模擬的運動模糊星圖，其他仿真條件見仿真條件設置.得到在不同的“種子”點誤差條件下，利用本文方法的星點提取誤差仿真結果如圖4、5所示.圖4中的星點提取成功率定義為ratio=Nsuc/Ntotal，其中Ntotal為用于仿真實驗的星圖總數，Nsuc為星圖上的星點被全部成功提取的星圖數量.

圖4 不同“種子”點誤差條件下的星點提取成功率

圖5 不同“種子”點誤差條件下err的均值

從本仿真實驗條件下，“種子”點誤差不超過3.5像素時，本文方法的星點提取成功率為100%，“種子”點誤差大于3.5像素時，星點提取成功率略有下降；仿真結果如圖5所示，星點提取誤差err的均值隨著“種子”點誤差的增大而略有增加，但增幅較小，上述仿真實驗結果說明本文方法對“種子”點誤差(即先驗信息誤差帶來的影響)具有較好的容錯能力.

3.2.2 仿真實驗2

根據t1時刻的光軸指向{(αi,δi)}得到100幅模擬的運動模糊星圖，其他仿真條件見仿真條件設置.利用本文方法直接對模糊星圖上星像質心進行提取和補償；作為比較，利用L-R方法和維納濾波方法對模糊星圖進行復原，然后對復原星圖利用塊掃描法[23]+質心法進行星點提取，分別簡記為Method-1和Method-2，L-R方法的迭代次數設置為10.得到仿真結果如圖6、7所示.

圖6 3種方法err值的對比

利用本文方法對100幅運動模糊星圖進行星點提取以及利用L-R方法以及維納濾波方法先對運動模糊星圖復原再進行星點提取，100幅星圖的星點提取精度如圖6所示，本文方法、Method-1以及Method-2這3種方法err的平均值分別為0.658、1.063、1.205像素，可以看出本文方法的星點提取誤差小于Method-1和Method-2；3種方法的耗時情況如圖7的所示，本文方法、Method-1和Method-2的平均耗時分別為0.063 2、0.555 3、0.112 0 s，可知本文方法的耗時小于Method-1和Method-2.從上述仿真實驗結果可知，相比于L-R方法和維納濾波方法，本文方法在快速性方面具有一定的優勢，星點提取誤差雖然比其他兩種方法有所降低，但誤差仍然較大.

3.2.3 仿真實驗3

為了更全面的分析星圖噪聲大小對星點提取的影響，研究本文方法的抗噪性能，設計如下仿真實驗.在模擬星圖上加入不同強度的高斯白噪聲，共設置10組不同強度的高斯白噪聲，高斯白噪聲均值為0，均方差分別為1,2,…,10(灰度值，星圖灰度值范圍為0～255)，其他仿真條件與仿真實驗2相同.得到仿真結果如圖8～10所示.

圖8 不同噪聲條件下3種方法的星點提取成功率對比

圖9 不同噪聲條件下3種方法耗時的對比

根據仿真結果圖8可知，隨著星圖噪聲增大，3種方法都受到了影響，Method-1的星點提取成功率下降幅度最大，Method-2次之，本文方法的星點提取成功率下降幅度最小，在仿真實驗中，星點提取成功率最低為0.8，可見本文方法具有較好的抗噪性能.從仿真結果圖9可知，3種方法的平均耗時隨著噪聲的變化，波動較小，相比于Method-1和Method-2，本文方法的速度更快.為了考察在不同噪聲條件下，當星圖上的星點被成功提取時3種方法的星點提取精度.在上述仿真實驗條件下，得到在不同噪聲條件下3種方法的星點提取誤差的仿真結果如圖10所示，圖10中的星點提取平均誤差為每組噪聲條件下err的均值.

圖10 不同噪聲條件下3種方法的星點提取誤差對比

從圖10的仿真結果可知，隨著噪聲的增大，3種方法的星點提取誤差均有逐漸增大的趨勢，但本文方法的星點誤差遠比其他兩種方法要小，以噪聲均方差等于10為例，本文方法的星點提取平均誤差為0.957像素，Method-1的星點提取平均誤差為1.784像素，Method-2的星點提取平均誤差為2.470像素.

根據仿真實驗2的結果，無論從星點提取成功率還是從星點提取平均誤差來討論，本文方法在不同噪聲水平下的表現都優于Method-1和Method-2，這說明本文方法具有較好的抗噪性能.

4 結 論

1)利用INS的輸出和星敏感器自身已經獲得的姿態先驗信息預測恒星在當前幀星圖上的粗略坐標，將其作為“種子”點，局部區域生長對“種子”點誤差具有較好的容錯能力，這在一定程度上減輕了先驗信息誤差對本文方法的影響.

2)建立了利用質心法提取運動模糊星圖上星點質心坐標的誤差修正模型.本文方法不需要對運動模糊星圖進行復原且不用對整幅星圖進行掃描，這大大節省了運算時間，有利于提高星敏感器數據更新率.

3)仿真實驗表明，本文方法的快速性和抗噪能力均優于常用的L-R方法和維納濾波方法，星點提取精度也得到了一定的提高.