陳 璐,王 昊,王本冬,金仲和
(浙江大學(xué)航空航天學(xué)院,杭州 310027)
近年來,微納衛(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)確定性能。
靜態(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.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。
當(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ù)。
靜態(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
地球紅外成像如圖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
式(6)~式(9)得到解畸變后的像點(diǎn)如圖12所示。圖中像點(diǎn)構(gòu)成的圓環(huán)即為截面圓環(huán)D,像點(diǎn)疏密不均,邊緣存在一些干擾點(diǎn)。
圖12 解畸變效果Fig.12 Distortion elimination effect
對(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)確定精度。
按照第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)確定效率。
針對(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à)值。