陳炳龍, 王 磊, 劉 幫, 周 衡

(中國科學(xué)院微小衛(wèi)星創(chuàng)新研究院, 上海 201306)

0 引 言

先進(jìn)天基太陽天文臺(tái)衛(wèi)星(advanced space-based solar observatory satellite, ASO-S)是我國第一顆綜合性太陽活動(dòng)觀測專用衛(wèi)星[1-2],科學(xué)目標(biāo)是“一磁兩爆”,即太陽磁場、太陽耀斑和日冕物質(zhì)拋射之間的關(guān)系,于2022年10月9日成功發(fā)射入軌。ASO-S衛(wèi)星搭載的科學(xué)載荷對(duì)衛(wèi)星平臺(tái)提出了高精度、高穩(wěn)定度對(duì)日指向控制要求。為此,衛(wèi)星總體要求姿軌控分系統(tǒng)的姿態(tài)確定精度在非對(duì)日方向(即衛(wèi)星-太陽連線垂直方向)優(yōu)于3″。同時(shí),為了保證長時(shí)間連續(xù)對(duì)日觀測,ASO-S衛(wèi)星軌道選擇了降交點(diǎn)地方時(shí)為6點(diǎn)的太陽同步軌道,建立如下對(duì)日觀測指向參考坐標(biāo)系:衛(wèi)星本體+Xb軸(即載荷太陽導(dǎo)行鏡(Sun guide telescope, SGT)的光軸方向)指向太陽中心、Yb軸平行于黃道平面。

基于星敏感器(star tracker, STR)和光纖陀螺(fiber optic gyro, FOG)的定姿算法是目前應(yīng)用最為廣泛的高精度姿態(tài)確定方法[3-5],因STR和FOG的測量精度較高而裝配在有高精度、高穩(wěn)定度需求的科學(xué)衛(wèi)星上。如2010年2月發(fā)射的太陽動(dòng)力學(xué)天文臺(tái)(solar dynamics observatory, SDO)衛(wèi)星裝備了2臺(tái)星敏感器和3臺(tái)兩軸角速度慣性測量設(shè)備[6]。經(jīng)過調(diào)研國內(nèi)衛(wèi)星使用的STR和FOG產(chǎn)品,ASO-S衛(wèi)星選用:3臺(tái)STR,姿態(tài)角測量精度優(yōu)于[3.5 3.5 28]″;2臺(tái)FOG,三軸角速度測量精度優(yōu)于5e-5°/s,常值零位漂移優(yōu)于0.2°/h。由于光學(xué)系統(tǒng)參數(shù)比較容易受到溫度影響,如星敏感器的主點(diǎn)位置、光纖陀螺的標(biāo)定因數(shù)等,都會(huì)影響敏感器的測量精度[7-9]。ASO-S衛(wèi)星將STR、FOG與載荷艙光學(xué)基準(zhǔn)板進(jìn)行隔熱安裝,利用SGT光學(xué)基準(zhǔn)立方鏡和帶有準(zhǔn)直功能的高精度經(jīng)緯儀確保STR和FOG裝配測量精度達(dá)到角秒[10-11]。同時(shí),熱控分系統(tǒng)對(duì)各單機(jī)進(jìn)行恒溫控制:STR安裝支架溫度控制在(0±3)℃,FOG和SGT安裝面溫度控制在(22±2)℃,使各單機(jī)獲得最佳性能。最后,為了減弱定姿使用STR切換產(chǎn)生的姿態(tài)控制偏差,在軌進(jìn)行3臺(tái)STR之間的安裝矩陣修正,并利用SGT數(shù)據(jù)修正STR安裝矩陣。通過注入指令修改相關(guān)姿控軟件參數(shù),將STR測量基準(zhǔn)統(tǒng)一至SGT測量坐標(biāo)系(即整星測量基準(zhǔn))。

盡管如此,基于STR和FOG的常規(guī)定姿算法仍難以滿足ASO-S的精度要求。為了獲得高精度和高魯棒性,國內(nèi)外學(xué)者對(duì)非線性濾波算法進(jìn)行了研究,如擴(kuò)展卡爾曼濾波器(extended Kalman filter, EKF)、無跡卡爾曼濾波器、粒子濾波器等[12-16]。但受限于星載計(jì)算機(jī)的處理能力,無法滿足非線性濾波算法的計(jì)算量。對(duì)于衛(wèi)星工程實(shí)踐來說,EKF定姿方法已被廣泛應(yīng)用[17-19]。

為實(shí)現(xiàn)高精度定姿,本文將載荷望遠(yuǎn)鏡中作為穩(wěn)像系統(tǒng)使用的SGT作為姿態(tài)敏感器使用,與STR、FOG測量數(shù)據(jù)一起組成濾波器測量系統(tǒng),設(shè)計(jì)姿態(tài)確定算法。通過數(shù)學(xué)仿真驗(yàn)證該方法的穩(wěn)定性和有效性。最后,使用ASO-S工程遙測數(shù)據(jù)驗(yàn)證該算法的姿態(tài)角估計(jì)精度,非對(duì)日方向定姿精度優(yōu)于0.2″。

1 系統(tǒng)方程

(1)

式中:ηstr為星敏感器測量誤差(白噪聲協(xié)方差矩陣為Qstr),記作ηstr～(0,Qstr);ω為衛(wèi)星本體系相對(duì)貫性系姿態(tài)角速度。

陀螺測量姿態(tài)角速度和常值零位漂移分別記作ωg和b(ωg=[ωgx,ωgy,ωgz]T,b=[bx,by,bz]T)。因此,衛(wèi)星慣性系角速度ω可以表示為

ω=ωg-b

(2)

(3)

(4)

2 測量方程

2.1 SGT測量模型

ASO-S科學(xué)載荷之一萊曼阿爾法太陽望遠(yuǎn)鏡使用高精度SGT作為穩(wěn)像系統(tǒng)。SGT利用兩組相對(duì)的光電探測器采集方位(zg)和俯仰方向(yg)的太陽邊緣信號(hào),測量原理如圖1所示。

圖1 SGT探測器焦平面示意圖Fig.1 Schematic diagram of SGT detector’s focal plane

圖1中A、B、C、D 4個(gè)光電二極管因光電效應(yīng)產(chǎn)生不同幅度的電流信號(hào),經(jīng)跨阻放大器后產(chǎn)生對(duì)應(yīng)的電壓信號(hào),分別為U1、U2、U3、U4,通過16位模數(shù)轉(zhuǎn)換器采集。太陽中心相對(duì)標(biāo)定的探測器焦平面中心偏移量與電壓信號(hào)滿足如下關(guān)系[20]:

(5)

根據(jù)偏移量計(jì)算出方位、俯仰方向(SGT探測器測量坐標(biāo)系與衛(wèi)星本體系三軸方向一致)偏差角:α=KyΔy、β=KzΔz(Ky和Kz為無量綱系數(shù))。

SGT為同軸光學(xué)系統(tǒng),太陽輻射首先經(jīng)過濾光片濾波后,進(jìn)入成像光學(xué)系統(tǒng),最后照射到二級(jí)管陣列探測器。影響SGT測量太陽中心偏差角的主要因素為光電轉(zhuǎn)換。通過太陽像輻照度分布函數(shù)I(θ)可獲得電流信號(hào),其中θ為日心角。當(dāng)太陽像中心偏移出光電二極管時(shí),其輸出特性呈現(xiàn)出明顯非線性[21-24]。因此,只有SGT在線性測量范圍時(shí),計(jì)算的太陽中心偏差角才能用于姿控分系統(tǒng)。此外,日心角隨著地球繞太陽公轉(zhuǎn)而變化,光電探測器上太陽像尺寸也會(huì)變化,需要在軌對(duì)無量綱系數(shù)進(jìn)行標(biāo)定。ASO-S在軌使用與SGT同光軸的白光望遠(yuǎn)鏡太陽圖像對(duì)SGT無量綱系數(shù)進(jìn)行標(biāo)定。入軌前由實(shí)驗(yàn)室以32′的太陽像進(jìn)行標(biāo)定,SGT的測量精度達(dá)到0.1″。

SGT測量原理可簡化成針孔成像模型[25-26],如圖2所示。

圖2 SGT測量原理圖Fig.2 Schematic diagram of SGT measurement principle

SGT焦距記作fgt,由此計(jì)算太陽矢量在SGT測量系的表達(dá)式為

(6)

太陽偏差角α、β與位置偏移量Δy、Δz滿足:Δy=f·tanα,Δz=f·tanβ。SGT線性測量范圍約為50″(tan 50″≈50″)。因此,Δy≈f·α、Δz≈f·β,rsun|gt改寫為

(7)

將慣性系太陽矢量Si轉(zhuǎn)換到衛(wèi)星本體系Sb:Sb=C(q)·Si(C(q)為四元數(shù)表示的姿態(tài)矩陣)。若Sb=[Sb_x,Sb_y,Sb_z]T、Si=[Si_x,Si_y,Si_z]T,則

式中:

Sb與rsun|gt平行,因此

可得到SGT測量的太陽中心偏差角α、β數(shù)學(xué)模型如下:

(8)

式中:ηgt為SGT的測量誤差(白噪聲協(xié)方差矩陣為Qgt),記作ηgt～ (0,Qgt)。

2.2 STR測量模型

STR測量方程可表示為

qstr=qis?qsb?ηstr

(9)

式中:?表示四元數(shù)乘法;qsb為STR的安裝四元數(shù);ηstr為測量誤差(白噪聲協(xié)方差矩陣為Qstr),記作ηstr～ (0,Qstr);qstr為計(jì)算出的本體系相對(duì)慣性系姿態(tài)四元數(shù)。

2.3 濾波器測量系統(tǒng)

選取STR測量值計(jì)算的衛(wèi)星本體系相對(duì)慣性系姿態(tài)四元數(shù)矢部qSTR_bv和SGT測量的偏差角α、β為測量量:y=[qSTR_b1,qSTR_b2,qSTR_b3,α,β]T,得到測量方程如下:

y=h(x)+v

(10)

3 EKF設(shè)計(jì)

(11)

(12)

式中:

濾波器輸入:由STR測量值計(jì)算的qSTR_bv,SGT測量值計(jì)算太陽中心偏差角α、β,FOG測量姿態(tài)角速度ωg,慣性系單位太陽矢量Si。

濾波器輸出:定姿四元數(shù)qnew,陀螺零漂估計(jì)bnew,定姿角速度ωnew。

按照如下步驟進(jìn)行迭代計(jì)算(Ts為計(jì)算步長):

步驟 1初始化

步驟 2均方誤差矩陣一步預(yù)測

步驟 3狀態(tài)一步預(yù)測

步驟 4濾波增益矩陣計(jì)算

Hk=H|qk,k-1

步驟 5狀態(tài)量更新

Xk=Xk,k-1+Kk(Zk-h|qk,k-1)

步驟 6估計(jì)均方誤差

Pk=(I-KkHk)Pk-1

步驟 7濾波器輸出

限幅:

qnew s

bnew_x/y/z>bmax→bnew_x/y/z=bmax

bnew_x/y/z<-bmax→bnew_x/y/z=-bmax

步驟 8定姿輸出

步驟 9更新濾波器參數(shù)

4 算法驗(yàn)證

4.1 定義姿態(tài)偏差角

(1) 定姿誤差

(2) 控制誤差

歐拉角形式:Δq轉(zhuǎn)換成歐拉角ΔΦ、ΔΘ、ΔΨ。

(3) 控制誤差與SGT指向誤差間的偏差角(在軌無法由真實(shí)姿態(tài)獲得定姿誤差,可用此偏差角反映定姿精度)

z軸偏差角:ΔΨ+α;

y軸偏差角:ΔΘ-β。

4.2 數(shù)學(xué)仿真條件

(1) 時(shí)間設(shè)置

開始時(shí)間:2022年5月10日1:40:00.000;

仿真時(shí)長:8 000 s;

仿真步長:0.125 s。

(2) 軌道參數(shù)

a=7.081 4×106m;e=9.010×10-4;i=1.716 5 rad;

Ω=2.402 2 rad;ω=4.278 7 rad;θ=3.891 5 rad。

(3) 初始姿態(tài)參數(shù)

姿態(tài)四元數(shù):q0=[0.891 8, 0.185 0,-0.083 8, 0.404 3]T;

姿態(tài)角速度:ω0=[0.0, 0.0, 0.0]T°/s;

轉(zhuǎn)動(dòng)慣量:Jb=diag([760.1, 537.7, 537.6])kg·m2。

(4) 敏感器參數(shù)

陀螺測量誤差:5×10-5°/s;

陀螺零位漂移:0.2°/h;

星敏測量誤差:[3.5″, 3.5″, 28″];

SGT測量誤差:0.1″;

SGT無量綱系數(shù):Ky=120.3、Kz=120.8;

陀螺安裝矩陣:CGYRO=diag([1,1,1]T);

星敏安裝矩陣(無偏差):

星敏安裝偏差(歐拉角):δESTR A=[30 15 30]″,δESTR B=[-15 30 -30]″,δESTR C=[30 -30 15]″;

星敏安裝矩陣(修正后):

(5) 執(zhí)行器參數(shù)

反作用輪最大力矩:0.215 N·m;

反作用輪最大摩擦力矩:0.025 N·m;

反作用輪類型:力矩輪;

反作用輪安裝方式:四斜裝“金字塔”構(gòu)型;

反作用輪力矩指令精度:12位數(shù)模轉(zhuǎn)換器;

磁力矩器最大磁矩:100 A·m2;

磁力矩器類型:開關(guān)式。

(6) 太陽模型

衛(wèi)星零時(shí)對(duì)應(yīng)的儒略日J(rèn)D0:

JD0=2 458 484.166 666 666 666 666 7

由衛(wèi)星時(shí)間積秒T計(jì)算JDUTC:

JDUTC=JD0+T/864 00

UTC時(shí)間2000年1月1日12:00:00對(duì)應(yīng)的儒略日:

JD2000=2 451 545.0

衛(wèi)星時(shí)間積秒相對(duì)2000年1月1日12時(shí)的儒略世紀(jì)數(shù)TRL:

TRL=(JDUTC-JD2000+69.184/864 00)/36 525

太陽軌道傾角:

Is=0.409 092 804 2-0.000 226 965 5TRL

太陽平近點(diǎn)角:

Ms=6.240 059 966 7+628.301 955 151 5TRL

太陽真黃經(jīng):

Us=4.895 062 993 9+628.307 584 536TRL+

0.033 416 073 9sinMs

J2000慣性系太陽位置單位矢量:

Si=[cosUs,sinUscosIs,sinUssinIs]T

考慮光行時(shí)和光行差效應(yīng)對(duì)觀測的影響:用迭代計(jì)算法對(duì)光行時(shí)進(jìn)行修正,用洛倫茲變換方法修正光行差,最后得到太陽位置單位矢量:Si_ls。

(7) 北黃極模型

由J2000慣性系轉(zhuǎn)換至任意JD歷元黃道坐標(biāo)系的轉(zhuǎn)換矩陣為

式中:Rz(*)表示繞任意坐標(biāo)系z(mì)軸轉(zhuǎn)動(dòng)某角度的姿態(tài)轉(zhuǎn)換矩陣;Rx(*)表示繞任意坐標(biāo)系x軸轉(zhuǎn)動(dòng)某角度的姿態(tài)轉(zhuǎn)換矩陣;ΨJ2000、φJ(rèn)2000、γJ2000為3個(gè)歲差角,單位為(″),表達(dá)式為

其中,TTT是自標(biāo)準(zhǔn)歷元J2000.0起算的(TT時(shí)刻)儒略世紀(jì)數(shù)(JDTT=JDUTC+69.184/864 00):

因此,J2000慣性系下北黃極單位矢量為

(8) 對(duì)日觀測導(dǎo)引律

定義衛(wèi)星對(duì)日觀測指向參考坐標(biāo)系的三軸在J2000慣性系如下:

(13)

因此,J2000慣性系轉(zhuǎn)換至對(duì)日觀測指向參考坐標(biāo)系的轉(zhuǎn)換矩陣可表示為

將姿態(tài)轉(zhuǎn)換矩陣計(jì)算成四元數(shù),即得到期望指向四元數(shù)qcmd。

(9) 比例-積分-微分(proportion integration defferentiation, PID)控制器參數(shù)

Kp=[24.22,17.85,17.61]T

Ki=[0.3,0.3,0.3]T

Kd=[130.10,95.88,94.58]T

(10) EKF參數(shù)

Q=diag([3.32e-13,3.32e-13,3.32e-13,

4.6e-15,4.6e-15,4.6e-15]T)

R=diag([5.88e-10,5.88e-10,5.88e-10,

2.35e-13,2.35e-13]T)

P0=diag([5.88e-10,5.88e-10,5.88e-10,

3.32e-13,3.32e-13,3.32e-13]T)

(11) 環(huán)境干擾力矩

重力梯度力矩Tg:

式中:μ為地球引力常數(shù),3.986 004 36×1014m3/s2;Z0為衛(wèi)星質(zhì)心相對(duì)地心矢徑的單位矢量;r為衛(wèi)星質(zhì)心至地心的距離。

剩磁力矩Tm:

Tm=Mm×B

式中:Mm為衛(wèi)星剩磁矩;B為衛(wèi)星所處位置的地球磁場向量。ASO-S衛(wèi)星對(duì)整星進(jìn)行了剩磁距補(bǔ)償,使整星剩磁矩<1 A·m2。

4.3 數(shù)學(xué)仿真結(jié)果

仿真共經(jīng)歷8 000 s (約1.3倍軌道周期),前1 000 s用于濾波器收斂。仿真結(jié)果如下:圖4為角速度估計(jì)誤差,圖5為姿態(tài)角估計(jì)誤差,圖6為導(dǎo)行鏡數(shù)據(jù)計(jì)算的指向偏差角(即衛(wèi)星實(shí)際指向誤差),圖7是由定姿結(jié)果計(jì)算的控制誤差與導(dǎo)行鏡計(jì)算指向誤差間的偏差角。

圖4 基于STR+FOG+SGT EKF的角速度估計(jì)值Fig.4 Estimated values of angular velocity by EKF basedon STR, FOG, and SGT

圖5 基于STR+FOG+SGT EKF的姿態(tài)角估計(jì)值Fig.5 Estimated values of attitude angle by EKF basedon STR, FOG, and SGT

圖6 由SGT測量值轉(zhuǎn)換的偏差角曲線Fig.6 Curve of deviation angles transformed frommeasurements of SGT

圖7 控制誤差與SGT計(jì)算指向誤差間的偏差角曲線Fig.7 Curve of deviation angles between controlerrors and pointing errors of SGT

以1 000 s至8 000 s數(shù)據(jù)計(jì)算角速度和姿態(tài)角估計(jì)精度,如表1和表2所示:基于SGT的EKF姿態(tài)角速度估計(jì)精度優(yōu)于1.5e-4°/s;x軸(對(duì)日指向)姿態(tài)角估計(jì)誤差優(yōu)于2.97″ (3σ),y、z軸(對(duì)日指向的垂直方向)姿態(tài)角估計(jì)誤差優(yōu)于0.12″ (3σ)。由SGT測量值計(jì)算太陽中心指向偏差精度:優(yōu)于1.0″(3σ)?？刂普`差與SGT指向誤差間的偏差角精度:優(yōu)于0.15″(3σ)。

表1 角速度估計(jì)精度Table 1 Estimated precision of angular velocity ((°)/s)

表2 姿態(tài)角估計(jì)精度Table 2 Estimated precision of attitude angle (″)

綜上所述,基于STR、FOG和SGT的EKF定姿系統(tǒng)可穩(wěn)定收斂并有效提高衛(wèi)星三軸姿態(tài)角估計(jì)精度:對(duì)日指向的估計(jì)精度優(yōu)于3.0″ (3σ);對(duì)日指向垂直方向的估計(jì)精度優(yōu)于0.12″ (3σ)。

4.4 在軌數(shù)據(jù)計(jì)算結(jié)果

使用ASO-S衛(wèi)星在軌工程遙測(UTC時(shí)間2023-03-06 09:08:30開始的3 350 s數(shù)據(jù)):姿控實(shí)時(shí)快速包的期望姿態(tài)四元數(shù)、慣性系姿態(tài)角速度、慣性系單位太陽矢量,星敏姿態(tài)數(shù)據(jù)包的3臺(tái)星敏感器數(shù)據(jù)有效性和星敏四元數(shù),導(dǎo)行鏡串口數(shù)據(jù)包的導(dǎo)行鏡y、z方向數(shù)據(jù),進(jìn)行算法驗(yàn)證。其中,導(dǎo)行鏡計(jì)算的指向偏差角如圖8所示。

圖8 在軌SGT測量值轉(zhuǎn)換的偏差角曲線Fig.8 Curve of deviation angles transformed fromon-orbit measurements of SGT

將ASO-S衛(wèi)星在軌工程數(shù)據(jù)輸入至本文設(shè)計(jì)的EKF至收斂,計(jì)算ΔΨ與SGT輸出指向偏差角α之和、ΔΘ與SGT輸出指向偏差角β之差,如圖9所示。由500 s至3 000 s間數(shù)據(jù)計(jì)算z、y軸偏差角誤差分別為0.161 5″ (3σ)和0.177 0″ (3σ)。使用星上進(jìn)行常規(guī)EKF濾波后的工程數(shù)據(jù)進(jìn)行相同計(jì)算,z、y軸偏差角如圖10所示,z、y軸偏差角誤差分別為4.400 7″ (3σ)和-11.287 1″ (3σ)。

圖9 以ASO-S衛(wèi)星在軌數(shù)據(jù)進(jìn)行EKF濾波后計(jì)算的z軸和y軸控制誤差與SGT計(jì)算指向誤差間的偏差角Fig.9 Deviation angles in axis z and y between controlerrors calculating by EKF and pointing errorsof SGT with ASO-S on-board data

圖10 ASO-S衛(wèi)星在軌進(jìn)行常規(guī)EKF濾波后的z軸和y軸控制誤差與SGT計(jì)算指向誤差間的偏差角Fig.10 Deviation angles in axis z and y between control errorscalculating by ASO-S on-board ordinaryEKF and pointing errors of SGT

綜上所述,以ASO-S衛(wèi)星在軌工程遙測輸入本文設(shè)計(jì)的濾波器,驗(yàn)證了該EKF算法可以快速收斂并穩(wěn)定輸;與常規(guī)EKF算法相比,該EKF算法顯著提高了對(duì)日指向垂直方向的定姿精度:由定姿結(jié)果計(jì)算指向誤差與SGT輸出真實(shí)指向誤差間的偏差角優(yōu)于0.2″(3σ)。

5 結(jié) 論

本文通過星敏、陀螺和SGT建立姿態(tài)測量系統(tǒng),以四元數(shù)表示的姿態(tài)運(yùn)動(dòng)學(xué)為系統(tǒng)方程,以四元數(shù)矢部和陀螺常值零位漂移為狀態(tài)變量,設(shè)計(jì)EKF定姿算法?？紤]敏感器安裝誤差,以比例-積分-微分控制器進(jìn)行數(shù)學(xué)仿真,結(jié)果表明:基于SGT的EKF姿態(tài)確定系統(tǒng)可穩(wěn)定收斂并有效提高衛(wèi)星三軸姿態(tài)角估計(jì)精度;對(duì)日指向優(yōu)于3.0″(3σ),對(duì)日指向垂直方向優(yōu)于0.15″(3σ)。最后,將ASO-S衛(wèi)星在軌工程遙測數(shù)據(jù)進(jìn)行基于SGT的EKF濾波,與常規(guī)EKF算法對(duì)比可知:本文設(shè)計(jì)的定姿算法可使對(duì)日指向垂直方向的姿態(tài)角估計(jì)精度優(yōu)于0.2″(3σ)。