褚理想， 樊巧云

(北京航空航天大學儀器科學與光電工程學院, 北京 100083)

姿態(tài)估計是衛(wèi)星上最重要的測量之一。目前常用的姿態(tài)傳感器有太陽敏感器、地磁、星敏感器、地球敏感器和GPS等，這些傳感器有其各自的優(yōu)點和使用限制。衛(wèi)星上常裝配多種敏感器，采用多種方案相互備份和補充，從而保證測量的可靠性和魯棒性。光電二極管[1]是一種成本低廉、體積更小的感光傳感器，根據(jù)其輸出電壓可以感測太陽矢量方向，可以作為低精度的太陽敏感器。

太陽敏感器與地球敏感器結合[2-3]，可以確定衛(wèi)星三軸姿態(tài)。然而基于光電二極管的太陽敏感器，其輸出的電壓是太陽光和地球反照光等雜光輻照的疊加，如果忽略地球反照光的影響，直接將測量電壓轉(zhuǎn)換為太陽矢量，會導致轉(zhuǎn)換結果存在較大誤差，進而影響姿態(tài)估計的精度。為此需要考慮地球反照光，建立精確的光電二極管量測模型，才能實現(xiàn)高精度的衛(wèi)星姿態(tài)估計。

文獻[4]提出了地球反照光的數(shù)學模型及其計算方法，該方法需要實時的球面積分運算，計算較為復雜。文獻[5-8]對其進行了簡化。文獻[5]假設地球表面的反照率僅與地球緯度有關，與地球經(jīng)度無關，建立地球反照光的查找表并將其擬合成曲面，在已知衛(wèi)星姿態(tài)和衛(wèi)星軌道位置的情況下可以快速計算地球反照光影響。這種方法前期計算較為復雜且僅適用于單一軌道。文獻[6]采用隨機實驗的方法得到地球反照光統(tǒng)計特性，然后根據(jù)統(tǒng)計的均值和方差補償太陽矢量模型，這種方法對太陽矢量的動態(tài)范圍有很大限制。文獻[7-8]假設地球反照光為與衛(wèi)星天底方向反向平行的單一矢量，簡化地球反照光的計算，但仍存在計算復雜和精度較低的問題。

本文結合地球反照光模型，提出光電二極管動態(tài)偏置混合高斯噪聲的量測模型，并應用滑窗估計和隨機權重策略動態(tài)更新地球反照光偏置和噪聲方差。在權重更新過程中，引入Huber影響函數(shù)處理量測殘差的異常值，從而提高了算法魯棒性和參數(shù)估計精度。同時，針對每個光電二極管的地球反照光分布不一，采用多比例因子分別估計每個光電二極管的地球反照光影響。基于本文建立的地球反照光校正的光電二極管太陽敏感器量測模型，結合地球敏感器和陀螺等傳感器，采用無跡卡爾曼濾波(Unscented Kalman Filter,UKF)算法[9]實現(xiàn)了衛(wèi)星的高精度姿態(tài)估計。

1 地球反照光校正的光電二極管量測模型

1.1 光電二極管工作原理

光電二極管的輸出電壓和太陽矢量與光電二極管敏感軸夾角的余弦成正比，根據(jù)電壓輸出，可以感測太陽方位。光電二極管常安裝于衛(wèi)星表面，易受周圍地球反照光等雜光的干擾。圖1為地球反照光幾何示意圖，照射在積分區(qū)域(太陽照射與衛(wèi)星可視方位區(qū)域交集)的太陽光，經(jīng)地球表面漫反射后，共同作用于衛(wèi)星方向。地球反照光的計算不僅與地球表面的反照率有關，還和太陽、衛(wèi)星、地球相對位置有關，此外，地球反照光還與衛(wèi)星本身的姿態(tài)有關。因此地球反照光的分布函數(shù)較為復雜，其計算較為繁瑣。

地球反照光的影響較大，可達到直射太陽光總量的20%～30%，不能簡單地忽略其影響?？紤]地球反照光的影響，單個光電二極管的輸出電壓V可以表示為[10]

V=Vd+Va+vV

(1)

(2)

圖1 地球反照光幾何示意圖Fig.1 Geometric sketch of the earth’s albedo

Va=

(3)

式中：Vd為太陽照射分量；Va為地球反照光分量；vV為模型誤差的零均值高斯白噪聲；n=[cosφcosθcosφsinθsinφ]T為光電二極管敏感軸的單位矢量方向，可由安裝高度角φ和方位角θ表示；sb是太陽矢量在本體坐標系下的表示，可由在慣性坐標下太陽矢量經(jīng)姿態(tài)矩陣變換得到；ψ為光電二極管半視場角；A為地球表面衛(wèi)星可視區(qū)域和太陽照射的交集區(qū)域；α為地球表面dA的反照率；nA為地球表面面元dA的法向單位矢量；s⊕為地球到太陽的矢量方向；rAB為地球表面面元dA到衛(wèi)星的矢量方向；B表示衛(wèi)星當前軌道位置；S表示在地球陰影區(qū)的軌道。

1.2 地球反照光校正

式(3)中地球反照光分量Va計算需要球面積分，運算復雜，難以在實際中得到應用。本文將地球反照光設成一動態(tài)偏置項，補償在光電二極管的量測模型中，同時將動態(tài)偏置估計的誤差和光電二極管本身量測的誤差統(tǒng)一為混合高斯模型，建立光電二極管的動態(tài)偏置混合高斯量測模型：

Vk=Vd,k+rk+vk

(4)

式中：下標k表示采樣時間序列；rk為地球反照光動態(tài)偏置；vk為非高斯噪聲，可用式(5)表示:

p(vk)=(1-ε)pN(vk)+εqN(vk)

(5)

其中：pN(vk)為已知的均值為0、方差為R1,k的高斯噪聲概率密度函數(shù)；qN(vk)為未知的均值為0、方差為R2,k的污染噪聲概率密度函數(shù)；參數(shù)ε∈(0,1)為污染系數(shù)，用來控制污染噪聲的強弱，本文取ε=0.05；p(vk)為整體量測噪聲vk的概率密度函數(shù)，其方差為Rk，滿足Rk=(1-ε)R1,k+εR2,k。

單個光電二極管僅能測量一個太陽矢量分量，至少需要有3個有效的光電二極管量測值才能求取完整的太陽矢量。對于多個光電二極管的量測模型，可以表示為

(6)

矢量表示形式為

Vk=Vd,k+rk+vk=Hsb,k+rk+vk

(7)

2 模型參數(shù)的在線估計和更新方法

式(7)中地球反照光校正的光電二極管量測模型的地球反照光偏置項rk和噪聲方差Rk是未知的，需要在衛(wèi)星姿態(tài)估計過程中在線估計并動態(tài)更新。本文采用滑窗估計和隨機權重策略動態(tài)更新模型參數(shù)，如圖2所示，假設模型參數(shù)在窗口采樣時間序列k-1,k-2,…,k-N范圍內(nèi)不變，對于歷史數(shù)據(jù)計算的模型參數(shù)賦予不同的權重，獲取當前時刻模型估計的參數(shù)。由式(7)，結合滑窗估計和隨機權重算法，此時可以得到動態(tài)偏置和噪聲方差的估計公式分別為

(8)

(9)

隨著量測噪聲的統(tǒng)計變化，量測殘差向量ek-j將存在偏置，它的幅值將會增加[11]，可以選用量測殘差的模值作為權重，如式(10)和式(11)所示：

wj∝|ek-j|j=1,2,…,N

(10)

(11)

然而光電二極管噪聲分布較為復雜，當出現(xiàn) 殘差較大的異常值時，直接采用殘差向量作為模值會引起權重分配不合理，進而導致偏差和方差估計產(chǎn)生較大誤差。基于Huber影響函數(shù)的魯棒技術可以有效地處理非高斯噪聲的情況，其更改了后驗噪聲方差矩陣，降低了對異常值的靈敏度。因此，本文引入Huber影響函數(shù)處理量測殘差的異常值，其表達式如下：

圖2 滑窗估計和隨機權重算法示意圖Fig.2 Schematic diagram of sliding window estimation and random weighting algorithm

(12)

(13)

式中：引入Tk使ηk滿足關于概率密度對稱和邊緣概率密度條件[12]。ψ(·)是Huber函數(shù)，其表達式為

(14)

其中:sgn(·)是符號函數(shù)；kε的取值和污染系數(shù)ε有關。

另外，考慮到每個光電二極管的地球反照光的分布情況不同，當所有的光電二極管都采用統(tǒng)一的權重系數(shù)會降低對地球反照光的跟蹤特性，因此本文采用多重比例因子分別估計地球反照光對每個光電二極管的影響，此時可得

(15)

(16)

式中：diag(·)表示將一個向量轉(zhuǎn)換為對角矩陣。

3 基于UKF的衛(wèi)星姿態(tài)估計

3.1 姿態(tài)估計狀態(tài)方程

衛(wèi)星姿態(tài)運動學[13]可以用四元數(shù)表示為

(17)

(18)

(19)

陀螺常用的數(shù)學模型為[15]

(20)

(21)

其中：δ(t-τ)為Dirac delta函數(shù)。

此時可以建立衛(wèi)星姿態(tài)估計的連續(xù)狀態(tài)方程為

(22)

δp=f[δqv/(a+δq4)]

(23)

式中：a為0到1區(qū)間的參數(shù)；f為放大因子。a=0 和f=1表示的是Gibbs向量，a=1和f=1表示的修正羅德里格斯向量。從δp到δq的逆變換為

(24)

δqv=f-1(a+δq4)δp

(25)

3.2 姿態(tài)估計量測方程

本文采用光電二極管和地球敏感器兩種姿態(tài)傳感器，需要將其測量值融入到量測方程中。光電二極管的量測方程另一種表達形式為

Vk=Hsb+rk+vk=HA(q)sref+rk+vk

(26)

式中：sref為太陽矢量在慣性下的表示，可以查找星歷表獲得。

地球敏感器測量天底方向，其量測模型可以表示為

bk=A(q)rearth+εk

(27)

結合2個傳感器測量模型可得量測方程：

(28)

3.3 UKF算法實現(xiàn)

當已知狀態(tài)更新的狀態(tài)方程模型，且建立了狀態(tài)和量測方程的噪聲和誤差統(tǒng)計模型，卡爾曼濾波方法采用遞推的方式，從量測信息中實時提取出被估計量信息并存儲在估計值中[17]。UKF是對線性卡爾曼濾波的改進，其不需要對狀態(tài)方程和量測方程線性化，常用于姿態(tài)估計等非線性濾波算法中。UKF通過sigma點捕獲系統(tǒng)真實的均值和方差，其精度可以達到泰勒展開式三階近似。

定義離散系統(tǒng)的非線性狀態(tài)方程和量測方程為

(29)

(30)

(31)

4 仿真校驗

4.1 仿真條件

n1=[cos(10°)cos(72°) cos(10°)sin(72°)

sin(10°)]T

n2=[cos(10°)cos(107°) cos(10°)sin(107°)

sin(10°)]T

n3=[cos(-20°)cos(90°) cos(-20°)sin(90°)

sin(-20°)]T

圖3為3個光電二極管的理想電壓和量測電壓的對比圖，地球反照光影響為量測電壓與理想電壓的差值。由圖中可以看出，對于同一個光電二極管，在軌道的不同時段,地球反照光與太陽直射光的比值不同,比值大約為0～25%。對于不同光電二極管，每個光電二極管的地球反照光分布情況不同。此外，在衛(wèi)星剛出背光面和剛?cè)氡彻饷鏁r，理想光電二極管和實際光電二極管的電壓相差不大，地球反照光較弱，可以選擇衛(wèi)星剛出背光面時刻作為太陽矢量估計或者姿態(tài)估計的初始時刻。

4.2 仿真結果

為了量化仿真結果，選取三軸歐拉角的模值作為評判指標：

(32)

式中：θk、φk和ψk分別為橫滾角、俯仰角和偏航角。為了保證結果的可靠性，參數(shù)使用蒙特卡羅仿真50次?；诒疚慕⒌牡厍蚍凑展庑Ｕ墓怆姸O管太陽敏感器量測模型，結合地球敏感器和陀螺等傳感器，采用UKF算法進行衛(wèi)星姿態(tài)估計。

圖4對比了固定權重(如均值)、量測殘差模值、量測殘差經(jīng)Huber影響函數(shù)處理后的模值、本文方法等4種權重選取策略的效果。從整體來看，在初始三軸衛(wèi)星姿態(tài)1.5°左右時，采用本文 建立的地球反照光校正模型和UKF算法，三軸姿態(tài)精度可以很快的收斂到0.5°甚至更高精度，驗證了本文簡化地球反照光模型具有一定的可行性。比較不同權重選擇策略，可以看出采用本文方法精度較高，三軸姿態(tài)精度可以達到0.2°～0.3°。

圖3 光電二極管1、2和3的理想電壓與量測電壓對比Fig.3 Comparison of ideal voltage and measured voltage for photodiode 1，2 and 3

圖4 衛(wèi)星三軸姿態(tài)估計誤差對比Fig.4 Comparison of satellite three-axis attitude estimation error

此外，還可以以地球反照光建模的動態(tài)偏置電壓估計精度作為評價標準，地球反照光估計精度越高，光電二極管測量的太陽矢量精度就越高，進而姿態(tài)的估計精度越高。單個光電二極管的偏置估計誤差可以通過多次蒙特卡羅方法求均值獲得。表1為光電二極管偏置估計誤差均方根(RMS)結果，可以看出，本文權重選取策略可以有效提高偏置估計精度。圖5為3個光電二極管的偏置估計誤差。從圖5(b)可以明顯看出,偏置估計的誤差隨時間推移而明顯減小，偏置估計的精度越來越高。

表1 光電二極管偏置估計誤差RMSTable 1 RMS of photodiode bias estimation error mV

圖5 光電二極管1、2和3的偏置估計誤差Fig.5 Bias estimation errors of photodiode 1,2 and 3

5 結 論

針對地球反照光的影響，本文建立了光電二極管動態(tài)偏置混合高斯量測模型，簡單有效，通過對地球反照光的準確估計，提高了光電二極管對太陽矢量的測量精度。實驗表明：

1) 應用地球反照光校正的光電二極管和地球敏感器組合定姿，可以消弱地球地球反照光的干擾，快速提高三軸姿態(tài)精度，精度可以達到0.2°～0.3°。

2) 在應用滑窗估計和隨機權重估計量測模型參數(shù)過程中，采用多比例因子和Huber影響函數(shù)的權重處理方法，可以有效提高地球反照光動態(tài)偏置電壓的估計精度。

本文提出的地球反照光校正方法，可以推廣應用于光電二極管與其他姿態(tài)傳感器(如地磁)的組合定姿。