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

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

0 引 言

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

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

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

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

1 算法分析

1.1 測(cè)量模型

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

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

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

y=fθ

(1)

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

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

1.2 成像點(diǎn)預(yù)處理

1.2.1邊緣腐蝕提取

由于地球紅外成像是一個(gè)圓面，后續(xù)計(jì)算是基于成像邊緣進(jìn)行的，故首先要進(jìn)行邊緣提取操作。這里采用Sobel算子強(qiáng)邊沿檢測(cè)和腐蝕結(jié)合的方式。具體而言，Sobel算子圖像領(lǐng)域?yàn)?/p>

(2)

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

計(jì)算梯度如下

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

(3)

(4)

二值圖像腐蝕即對(duì)圖像A采用向量減法運(yùn)算

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

(5)

式中：B(X,Y)為位于(X，Y)的結(jié)構(gòu)元素像素區(qū)域，運(yùn)算AΘB表示:與A相應(yīng)位置的像素值一致時(shí)，A相應(yīng)位置的像素點(diǎn)保留[15]。

1.2.2解畸變

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

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

(6)

圖2 像素點(diǎn)解畸變Fig.2 Elimination of pixel distortion

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

(7)

(8)

z=Rcos(θ)

(9)

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

1.3 目標(biāo)角計(jì)算

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

圖3 算法流程圖對(duì)比Fig.3 Comparison of algorithm flow charts

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

1.3.1間隔向量積

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

圖4 像素點(diǎn)映射至空間位置Fig.4 Pixel mapping spatial position

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

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

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

從這三個(gè)集合中選取三個(gè)點(diǎn)ak,bk和ck，得到兩個(gè)向量Vbak和Vbck，稱為空間向量對(duì)。順序取點(diǎn)，至少可以構(gòu)成n個(gè)空間矢量對(duì)，它們?cè)谕豢臻g圓環(huán)D上。由于向量Vbak和向量Vbck共面不共線，根據(jù)向量理論[17]，得到

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

(10)

式中：Vβk為向量積結(jié)果，Vβ是由n個(gè)Vβk向量構(gòu)成的集合。

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

(11)

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

1.3.2t Location-Scale分布擬合

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

圖6 5°時(shí)Bt內(nèi)樣本概率密度分布Fig.6 Probability density distribution of samples in Bt at 5°

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

(12)

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

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

(13)

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

(14)

(15)

(16)

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

(17)

(18)

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

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

(19)

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

(20)

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

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

(21)

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

1.3.3置信區(qū)間均值

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

(22)

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

2 試驗(yàn)測(cè)試

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

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

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

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

2)PC機(jī)發(fā)出角度控制信號(hào)，控制轉(zhuǎn)臺(tái)轉(zhuǎn)動(dòng)。

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

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

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

圖8 測(cè)試場(chǎng)景Fig.8 Test scenario

3 測(cè)試結(jié)果分析

3.1 邊緣提取

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

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

圖10 Sobel邊緣檢測(cè)Fig.10 Sobel edge detection

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

3.2 解畸變

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

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

3.3 間隔向量積與分布結(jié)果

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

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

之后通過t Location-Scale擬合，得到Bt樣本概率密度分布和擬合結(jié)果。例如，當(dāng)目標(biāo)β值為60°，30°和5°時(shí)，結(jié)果如圖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°以內(nèi)，Bt樣本分布左偏。t Location-Scale分布估計(jì)的參數(shù)，如表1所示。由表1可知，隨著β值的減小，其均值附近的概率密度顯著增加，分布趨向于集中。之后計(jì)算在1σ區(qū)間內(nèi)Bt樣本均值，如圖14所示。對(duì)每一張紅外成像做如上處理，可以得到目標(biāo)角度β的估計(jì)值。

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

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

3.4 算法結(jié)果比較

按照第2.2節(jié)所述的試驗(yàn)方法進(jìn)行測(cè)試。將新算法與常用的霍夫變換+最小二乘法[10]的實(shí)測(cè)結(jié)果比較，如圖16和圖17所示。可知當(dāng)測(cè)量角度大于40°時(shí)，新算法能夠得到平滑的測(cè)試結(jié)果。

圖16 原算法目標(biāo)角度估計(jì)結(jié)果Fig.16 Target angle estimation results of original algorithm

圖17 新算法目標(biāo)角度估計(jì)結(jié)果Fig.17 Target angle estimation results of new algorithm

圖18 ±10°內(nèi)原算法測(cè)量誤差Fig.18 Measuring error of original algorithm in ± 10°

圖19 ±10°內(nèi)新算法測(cè)量誤差Fig.19 Measuring error of original algorithm in ± 10°

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

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

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

4 結(jié) 論

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