基于太阳导行镜测量的高精度姿态确定算法

陈炳龙, 王 磊, 刘 帮, 周 衡

(中国科学院微小卫星创新研究院, 上海 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σ)。