陳 璐，王 昊，王本冬，金仲和

(浙江大學航空航天學院，杭州 310027)

0 引 言

近年來，微納衛(wèi)星技術得到了迅猛發(fā)展。它有著質量輕、體積小、功耗小的特點[1-2]，在航天應用領域極具潛力。目前，微納衛(wèi)星所執(zhí)行的任務種類和難度在不斷增加，例如要求微納衛(wèi)星實現(xiàn)大角度快速機動、編隊飛行、雙目視覺等[3-5]。而提升微納衛(wèi)星的機動能力是實現(xiàn)高難度任務的基礎。

微納衛(wèi)星機動能力取決于姿態(tài)敏感器和控制器的性能。姿態(tài)敏感器包括地球敏感器、太陽敏感器、星敏感器和陀螺儀等[6-7]。其中靜態(tài)紅外地球敏感器是微納衛(wèi)星較為常用且可靠的定姿敏感器[8-12]。

靜態(tài)紅外地球敏感器通過對地球凝視拍攝，計算成像矢量的方向，經(jīng)過坐標轉換，最終得到衛(wèi)星的對地姿態(tài)[9]。它相較于動態(tài)掃描式紅外地球敏感器而言，具有速度快、功耗小、壽命長和體積小等優(yōu)勢。文獻[10]中靜態(tài)紅外地球敏感器的姿態(tài)確定算法是將最小二乘法與霍夫變換結合，通過計算地球圓盤圓心來確定姿態(tài)角，姿態(tài)確定精度為0.1°。文獻[11] 也利用最小二乘迭代的方式求解非整圓弧的圓心。文獻[12]考慮了地球扁率和紅外輻射曲線隨季節(jié)變化的特點，利用補償和濾波的方式獲得地球切線邊緣。但以上算法都基于地球圓盤圓心的計算，這會占用較多計算資源，同時地球圓盤圓心坐標的微小偏差也會使姿態(tài)角計算產生明顯偏差。

本文提出了一種基于靜態(tài)全景紅外地球敏感器的微納衛(wèi)星姿態(tài)角確定算法。新算法采用間隔向量積與t Location-Scale分布置信區(qū)間均值結合的算法。相較于傳統(tǒng)的方法，新算法直接計算地球圓環(huán)的法向量，避免了計算圓環(huán)圓心位置的中間步驟，使得定姿精度提高50%，定姿速度提高87%，提升了微納衛(wèi)星的姿態(tài)確定性能。

1 算法分析

1.1 測量模型

靜態(tài)紅外地球敏感器的測量目標如圖1所示，地心矢量與成像軸Z軸的夾角β為測定的姿態(tài)角。

圖1 靜態(tài)紅外地球敏感器的成像示意圖Fig.1 Schematic diagram of static infrared earth sensor imaging

敏感器采用全景式魚眼紅外鏡頭，其成像規(guī)律為[13]：

y=fθ

(1)

式中：y為入射光線所成像素點與成像中心的距離，f為焦距，θ為入射光線與鏡頭法線夾角。

紅外敏感器的成像與θ正相關，故無論敏感器轉到什么方向，地球成像面均為E面。將E面投射到光學系統(tǒng)成像面，其坐標為(X，Y)，系統(tǒng)算法的目的是通過坐標(X，Y)得到β值。

1.2 成像點預處理

1.2.1邊緣腐蝕提取

由于地球紅外成像是一個圓面，后續(xù)計算是基于成像邊緣進行的，故首先要進行邊緣提取操作。這里采用Sobel算子強邊沿檢測和腐蝕結合的方式。具體而言，Sobel算子圖像領域為

(2)

式中：z1,z2,z3,…,z9為該位置像素點的灰度值。

計算梯度如下

[(z3+2z6+z9)-(z1+2z4+z7)]2}1/2

(3)

(4)

二值圖像腐蝕即對圖像A采用向量減法運算

AΘB={(X,Y)|B(X,Y)?A}

(5)

式中：B(X,Y)為位于(X，Y)的結構元素像素區(qū)域，運算AΘB表示:與A相應位置的像素值一致時，A相應位置的像素點保留[15]。

1.2.2解畸變

基于全景魚眼紅外鏡頭的成像特性[13]，將成像點(X，Y)映射至半徑為R的球面上，得到的邊沿點(x,y,z)構成一個截面圓環(huán)D，如圖2所示。

根據(jù)鏡頭成像公式(1)有

(6)

圖2 像素點解畸變Fig.2 Elimination of pixel distortion

根據(jù)物理成像關系，可以得到成像位置坐標與紅外成像三維方向之間的關系:

(7)

(8)

z=Rcos(θ)

(9)

式中：R為映射球面半徑，點(x,y,z)為光學成像點(X，Y)到球面的映射，它們構成一個空間截面圓環(huán)D。

1.3 目標角計算

當?shù)玫匠上駡A環(huán)后，通常利用圓心位置計算目標角度，方法包括：三點均值法，霍夫變換+最小二乘法和面積積分法[16]。本文提出間隔向量積與t Location-Scale分布置信區(qū)間均值結合的方法，與原有霍夫變換+最小二乘算法[10]對比，如圖3所示。

圖3 算法流程圖對比Fig.3 Comparison of algorithm flow charts

間隔向量積和t Location-Scale分布擬合是新算法與傳統(tǒng)算法的本質差別。一方面，間隔向量積采用直接計算圓環(huán)法向量的方法，避免了計算圓環(huán)圓心的中間步驟，也避免了霍夫變換對切半徑的估計和全局計算。減少了關鍵信息的損耗，節(jié)省計算資源。另一方面，t Location-Scale分布擬合可用于總體方差未知的情況，能較好地排除異常值，模擬尾部較重的數(shù)據(jù)分布和集中趨勢。由t Location-Scale分布擬合濾去干擾點，可以獲得準確的姿態(tài)角。

1.3.1間隔向量積

經(jīng)過解畸變處理后，光學系統(tǒng)成像平面上的像素點被投影到空間球面上，構成截面圓環(huán)D,如圖4所示。圓環(huán)D上每個像素點有一個三維坐標(x,y,z)。然而由于邊緣提取的誤差問題，有些像素點可能不在截面圓上，它們被稱為干擾點。

圖4 像素點映射至空間位置Fig.4 Pixel mapping spatial position

任意選取截面圓環(huán)D上三點a，b和c，可以構成兩個向量Vba和Vbc。由向量定理可知[15]，兩向量Vba和Vbc矢積得到的向量Vβ與原兩向量所在平面D垂直。則有Vβ與r平行，顯然Vβ與Z軸所成夾角即為目標值β。

具體而言，若將a,b,c三點選取的距離較近，其構成的兩向量將趨于共線，矢積所得垂向量趨于0，不便于后續(xù)計算，所以應間隔選取像點。將空間圓環(huán)D分為A段、B段和C段，表示為三個集合A={a1,a2,…,a3},B={b1,b2,…,b3}和C={c1,c2,…,c3}。各部分包含n個像點，如圖5所示。

圖5 分割截面圓環(huán)D上像點Fig.5 Partition of image points on ring D

從這三個集合中選取三個點ak,bk和ck，得到兩個向量Vbak和Vbck，稱為空間向量對。順序取點，至少可以構成n個空間矢量對，它們在同一空間圓環(huán)D上。由于向量Vbak和向量Vbck共面不共線，根據(jù)向量理論[17]，得到

Vβ={Vβk|Vβk=Vabk×Vack,k=1,2,3,…,n}

(10)

式中：Vβk為向量積結果，Vβ是由n個Vβk向量構成的集合。

Vβ中向量Vβk均垂直于圓環(huán)D，故計算β角有

(11)

式中：Bt為單張紅外成像得到的β角集合。

1.3.2t Location-Scale分布擬合

由于干擾點的存在，有些像點并不分布在同一平面D上，導致了Bt集合當中角度βk的差異。其典型的概率密度分布如圖6所示。

圖6 5°時Bt內樣本概率密度分布Fig.6 Probability density distribution of samples in Bt at 5°

由圖6可知，βk的分布存在偏度，即存在較重尾部。t Location-Scale分布可用于總體方差未知的情況，能夠模擬出尾部較重的數(shù)據(jù)分布，擬合數(shù)據(jù)的集中趨勢，同時排除異常值的干擾[18]。相比之下，受重尾的影響正態(tài)分布的均值容易偏離大多數(shù)樣本所在的區(qū)域。所以選擇t Location-Scale分布進行擬合，其概率密度函數(shù)如下

(12)

式中：Γ(·)表示Gamma方程，∞ <μ<+∞為位置參數(shù)，σ>0為尺度參數(shù)，ν>0為形狀參數(shù)。

t Location-Scale分布的參數(shù)可以通過極大似然函數(shù)估計得到[19]則其極大似然函數(shù)有

(13)

對(μ,σ,ν) 求μ和σ的偏導數(shù)，令其等于0有

(14)

(15)

(16)

式(14)至式(16)表明，參數(shù)間存在耦合，但仍可以通過迭代收斂的方法得到估計值，其迭代式為

(17)

(18)

式中：(i)表示第i次迭代。

迭代終止條件為前后兩次結果差值小于10-4。由于自由度ν未知，故進行極大似然估計的迭代

(19)

式中：φ(x)=dln(Γ(x))/dx,-φ(·)+ln(·)近似為[20]

(20)

由式(17)～式(20)可得t Location-Scale的μ和σ。

在具體工程應用中，實時通過t Location-Scale分布擬合，需占用較多計算資源，故簡化計算。由于Bt∈(0,90)°，故將區(qū)間(0,90)劃分為h個子區(qū)間(90(i-1)/h, 90i/h)，其中i=1,2,3,…,h，則估算μ有

(21)

地面試驗可測得的t Location-Scale分布在各個角度區(qū)間的方差，由此建立方差經(jīng)驗表，如表1所示，實際工程中可以通過μ查表得到σ。由于t Location-Scale分布的方差結果主要受紅外鏡頭固有性能的影響，在重復測試后，固定角度的方差估計結果沒有大的波動。所以方差經(jīng)驗表的引入在±10°內對衛(wèi)星姿態(tài)確定精度沒有明顯的影響，同時能減少迭代估計的運算量，從而加快運算速度。

1.3.3置信區(qū)間均值

由于Bt樣本較易出現(xiàn)異常值，直接利用整體均值μ估計β值，會受到干擾。故選取一定置信區(qū)間內的樣本均值進行估計。雖然置信區(qū)間較寬時，其置信水平較高，但是包含的干擾點較多。根據(jù)3σ準則可知，數(shù)據(jù)分布在(μ-σ,μ+σ)的概率為68.27%，是較為適中的接受域。最終估計β值為

(22)

式中：m為落入接受域中βk的個數(shù)。

2 試驗測試

靜態(tài)紅外地球敏感器的測試系統(tǒng)，如圖7所示。

圖7 靜態(tài)紅外地球敏感器測試系統(tǒng)Fig.7 Static infrared earth sensor test system

圖7中紅外地球敏感器測試系統(tǒng)由5部分組成：紅外地球模擬器、紅外地球敏感器、單軸高精度轉臺、DSP和PC機。測試方案如下：

1)紅外地球模擬器的直徑為38 cm，鏡頭距離模擬地球的距離為3 cm，模擬950 km軌道高度。

2)PC機發(fā)出角度控制信號，控制轉臺轉動。

3)將紅外地球敏感器置于高精度轉臺上，轉臺帶動紅外地球敏感器轉動，多角度拍攝地球模擬器。在(-10°，10°)范圍內，按照 0.1°為間隔拍攝，在(-50°，-10°)和(10°，50°)范圍，按照1°為間隔拍攝。得到不同角度下的地球成像。

4)DSP將成像數(shù)據(jù)返回給PC，PC機利用圖像處理算法分析數(shù)據(jù)，得出地球矢量信息。

根據(jù)測試系統(tǒng)開展測試，測試場景如圖8所示。

圖8 測試場景Fig.8 Test scenario

3 測試結果分析

3.1 邊緣提取

地球紅外成像如圖9所示，需將圖像中間60°盲區(qū)去除。之后由Sobel邊緣檢測得到強邊緣圖像如圖10所示。最后經(jīng)過二值圖像腐蝕，如圖11所示。

圖9 地球紅外成像Fig.9 Earth infrared imaging

圖10 Sobel邊緣檢測Fig.10 Sobel edge detection

圖11 二值圖像腐蝕Fig.11 Binary image corrosion

3.2 解畸變

式(6)～式(9)得到解畸變后的像點如圖12所示。圖中像點構成的圓環(huán)即為截面圓環(huán)D，像點疏密不均，邊緣存在一些干擾點。

圖12 解畸變效果Fig.12 Distortion elimination effect

3.3 間隔向量積與分布結果

對截面圓環(huán)D上的像點進行間隔向量積運算，以10°附近為例，可以得到運算結果如圖13所示。

圖13 10度的Bt樣本Fig.13 Bt samples at 10°

之后通過t Location-Scale擬合，得到Bt樣本概率密度分布和擬合結果。例如，當目標β值為60°，30°和5°時，結果如圖14、圖15和圖6所示。

圖14 60°Bt樣本概率密度分布Fig.14 Probability density distribution of samples in Bt at 60°

從圖14和圖15可以看出，Bt樣本分布呈現(xiàn)單峰特性，存在偏度，在60°附近，Bt樣本分布右偏；在30°以內，Bt樣本分布左偏。t Location-Scale分布估計的參數(shù)，如表1所示。由表1可知，隨著β值的減小，其均值附近的概率密度顯著增加，分布趨向于集中。之后計算在1σ區(qū)間內Bt樣本均值，如圖14所示。對每一張紅外成像做如上處理，可以得到目標角度β的估計值。

表1 t Location-Scale分布參數(shù)估計Table 1 t Location-Scale distribution parameter estimation

值得注意的是，以下因素會對新算法造成影響。首先，t Location-Scale分布的形狀參數(shù)過高，將使得分布接近于正態(tài)分布，減弱剔除干擾點的能力；反之分布曲線雙側尾部會過高，將減少算法的有效數(shù)據(jù)量。尺度參數(shù)過高會使計算包含太多干擾點；反之會減弱算法對位置參數(shù)不準確的容忍能力。其次，間隔向量積中像點距離不能隨機選取，否則將可能使像點構成的兩向量趨于共線，從而無法計算姿態(tài)角。另外，紅外地球成像的輪廓不能過腐蝕，否則將損耗有效像點，影響算法的姿態(tài)確定精度。

3.4 算法結果比較

按照第2.2節(jié)所述的試驗方法進行測試。將新算法與常用的霍夫變換+最小二乘法[10]的實測結果比較，如圖16和圖17所示?？芍敎y量角度大于40°時，新算法能夠得到平滑的測試結果。

圖16 原算法目標角度估計結果Fig.16 Target angle estimation results of original algorithm

圖17 新算法目標角度估計結果Fig.17 Target angle estimation results of new algorithm

圖18 ±10°內原算法測量誤差Fig.18 Measuring error of original algorithm in ± 10°

圖19 ±10°內新算法測量誤差Fig.19 Measuring error of original algorithm in ± 10°

由圖18和圖19可知，當目標角度在±10°以內時，原算法的誤差角度的標準差為0.06°，新算法為0.03°，測量精度提升了50%。進一步比較各區(qū)間內算法精度及計算速度如表2所示。

表2 在不同區(qū)間內算法的精度比較Table 2 Accuracy comparison of algorithms in different regions

由表2可知，新算法在各區(qū)間的測量精度均有所提高。另外，表2中列出了新算法工程版的性能，其單張紅外成像的平均處理時間相對于原算法縮短了87%，提升了靜態(tài)紅外地球敏感器姿態(tài)確定效率。

4 結 論

針對靜態(tài)紅外地球敏感器的姿態(tài)角度計算問題，新算法采用了間隔向量積與t Location-Scale置信區(qū)間均值結合的方法。相比于原算法中霍夫變換與最小二乘法結合的方法，新算法無需給定相對投影距離，減少了計算的不確定性。試驗表明，在±10°的視場內，新算法的姿態(tài)確定精度為0.03°，與原算法相比提高了50%，同時其計算時間縮短了87%。本研究提升了靜態(tài)紅外地球敏感器的測量精度和效率，對于微納衛(wèi)星的姿態(tài)確定有一定的應用價值。