基于太陽導行鏡測量的高精度姿態確定算法

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

(中國科學院微小衛星創新研究院, 上海 201306)

0 引 言

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

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

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

為實現高精度定姿,本文將載荷望遠鏡中作為穩像系統使用的SGT作為姿態敏感器使用,與STR、FOG測量數據一起組成濾波器測量系統,設計姿態確定算法。通過數學仿真驗證該方法的穩定性和有效性。最后,使用ASO-S工程遙測數據驗證該算法的姿態角估計精度,非對日方向定姿精度優于0.2″。

1 系統方程

(1)

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

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

ω=ωg-b

(2)

(3)

(4)

2 測量方程

2.1 SGT測量模型

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

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

圖1中A、B、C、D 4個光電二極管因光電效應產生不同幅度的電流信號,經跨阻放大器后產生對應的電壓信號,分別為U1、U2、U3、U4,通過16位模數轉換器采集。太陽中心相對標定的探測器焦平面中心偏移量與電壓信號滿足如下關系[20]:

(5)

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

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

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

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

SGT焦距記作fgt,由此計算太陽矢量在SGT測量系的表達式為

(6)

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

(7)

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

式中:

Sb與rsun|gt平行,因此

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

(8)

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

2.2 STR測量模型

STR測量方程可表示為

qstr=qis?qsb?ηstr

(9)

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

2.3 濾波器測量系統

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

y=h(x)+v

(10)

3 EKF設計

(11)

(12)

式中:

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

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

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

步驟 1初始化

步驟 2均方誤差矩陣一步預測

步驟 3狀態一步預測

步驟 4濾波增益矩陣計算

Hk=H|qk,k-1

步驟 5狀態量更新

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

步驟 6估計均方誤差

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更新濾波器參數

4 算法驗證

4.1 定義姿態偏差角

(1) 定姿誤差

(2) 控制誤差

歐拉角形式:Δq轉換成歐拉角ΔΦ、ΔΘ、ΔΨ。

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

z軸偏差角:ΔΨ+α;

y軸偏差角:ΔΘ-β。

4.2 數學仿真條件

(1) 時間設置

開始時間:2022年5月10日1:40:00.000;

仿真時長:8 000 s;

仿真步長:0.125 s。

(2) 軌道參數

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) 初始姿態參數

姿態四元數:q0=[0.891 8, 0.185 0,-0.083 8, 0.404 3]T;

姿態角速度:ω0=[0.0, 0.0, 0.0]T°/s;

轉動慣量:Jb=diag([760.1, 537.7, 537.6])kg·m2。

(4) 敏感器參數

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

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

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

SGT測量誤差:0.1″;

SGT無量綱系數: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) 執行器參數

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

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

反作用輪類型:力矩輪;

反作用輪安裝方式:四斜裝“金字塔”構型;

反作用輪力矩指令精度:12位數模轉換器;

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

磁力矩器類型:開關式。

(6) 太陽模型

衛星零時對應的儒略日JD0:

JD0=2 458 484.166 666 666 666 666 7

由衛星時間積秒T計算JDUTC:

JDUTC=JD0+T/864 00

UTC時間2000年1月1日12:00:00對應的儒略日:

JD2000=2 451 545.0

衛星時間積秒相對2000年1月1日12時的儒略世紀數TRL:

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

太陽軌道傾角:

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

太陽平近點角:

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

太陽真黃經:

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

0.033 416 073 9sinMs

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

Si=[cosUs,sinUscosIs,sinUssinIs]T

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

(7) 北黃極模型

由J2000慣性系轉換至任意JD歷元黃道坐標系的轉換矩陣為

式中:Rz(*)表示繞任意坐標系z軸轉動某角度的姿態轉換矩陣;Rx(*)表示繞任意坐標系x軸轉動某角度的姿態轉換矩陣;ΨJ2000、φJ2000、γJ2000為3個歲差角,單位為(″),表達式為

其中,TTT是自標準歷元J2000.0起算的(TT時刻)儒略世紀數(JDTT=JDUTC+69.184/864 00):

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

(8) 對日觀測導引律

定義衛星對日觀測指向參考坐標系的三軸在J2000慣性系如下:

(13)

因此,J2000慣性系轉換至對日觀測指向參考坐標系的轉換矩陣可表示為

將姿態轉換矩陣計算成四元數,即得到期望指向四元數qcmd。

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

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參數

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) 環境干擾力矩

重力梯度力矩Tg:

式中:μ為地球引力常數,3.986 004 36×1014m3/s2;Z0為衛星質心相對地心矢徑的單位矢量;r為衛星質心至地心的距離。

剩磁力矩Tm:

Tm=Mm×B

式中:Mm為衛星剩磁矩;B為衛星所處位置的地球磁場向量。ASO-S衛星對整星進行了剩磁距補償,使整星剩磁矩<1 A·m2。

4.3 數學仿真結果

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

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

圖5 基于STR+FOG+SGT EKF的姿態角估計值Fig.5 Estimated values of attitude angle by EKF basedon STR, FOG, and SGT

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

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

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

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

表2 姿態角估計精度Table 2 Estimated precision of attitude angle (″)

綜上所述,基于STR、FOG和SGT的EKF定姿系統可穩定收斂并有效提高衛星三軸姿態角估計精度:對日指向的估計精度優于3.0″ (3σ);對日指向垂直方向的估計精度優于0.12″ (3σ)。

4.4 在軌數據計算結果

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

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

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

圖9 以ASO-S衛星在軌數據進行EKF濾波后計算的z軸和y軸控制誤差與SGT計算指向誤差間的偏差角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衛星在軌進行常規EKF濾波后的z軸和y軸控制誤差與SGT計算指向誤差間的偏差角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衛星在軌工程遙測輸入本文設計的濾波器,驗證了該EKF算法可以快速收斂并穩定輸;與常規EKF算法相比,該EKF算法顯著提高了對日指向垂直方向的定姿精度:由定姿結果計算指向誤差與SGT輸出真實指向誤差間的偏差角優于0.2″(3σ)。

5 結 論

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