卫星姿态确定与三轴磁强计校正_黄琳

第29卷第3期2008年5月　宇　航　学　报Journal of AstronauticsVol .29May No .32008卫星姿态确定与三轴磁强计校正黄　琳,荆武兴(哈尔滨工业大学航天工程系,哈尔滨150001) 摘　要:讨论了三轴磁强计观测存在系统误差情况下的地磁姿态确定问题并提出几个可行算法。

首先,通过将三轴磁强计系统误差参数扩展到四元数和陀螺常漂组成的原状态量中,提出一个姿态确定扩维UKF 滤波器,接着利用姿态确定与参数辨识双重滤波器思想提出2个与姿态信息有关 无关的姿态确定双重UKF 滤波器。

最后,利用计算机仿真分析了扩维状态系统的可观性并检验了新滤波器的有效性与可行性。

关键词:地磁姿态确定;三轴磁强计;计算机仿真;近地卫星中图分类号:V448.2 文献标识码:A 文章编号:1000-1328(2008)03-0854-06收稿日期:2007-07-09;　修回日期:2007-12-030　引言基于三轴磁强计(Three -Axis Ma gnetometer ,简记为TAM )观测的近地卫星姿态确定方法[1-5,11]已经发展起来。

但是,大多数方法[1-4,11]都采用简单的地磁观测模型,仅仅考虑了随机观测噪声的影响,而没有完整考虑其它诸多系统误差因素的影响。

比较典型的系统误差影响因素有TAM 本身三个磁通线圈轴的非正交性、刻度因子调试误差、零点常偏,以及剩磁场干扰[5]和电磁设备工作时产生的附加干扰磁场[6]。

一些实际案例[7]也证实对TAM 进行在线校正以确保获得一定精度的姿态确定结果来说是非常必要的。

Crassidis 等人[8]提出三个适合在线校正T AM 系统误差(3个非正交因素、3个刻度因子误差和3个常偏置量)的参数滤波算法,即序贯中心化方法、普通EKF 和简单UKF 校正算法。

事实上,后两者的实质区别在于新息协方差阵的计算有所不同,但UKF 滤波器通常能够取得相对其它两个滤波器来说更高一些的估计精度和更快的收敛速度,而且对有色观测噪声也显示出良好的鲁棒性。

需要指出的是,这三个校正算法仅仅利用了地磁观测量的大小,而与其指向(或卫星姿态)无关,属于所谓的标量校正算法。

目前,利用非线性滤波理论解决姿态确定与TAM 在线校正这一组合估计问题还是一个新鲜课题。

本文作者曾经讨论了地磁姿态确定与地磁观测常偏置量校正两者组合的估计问题[5],该偏置量实际上是稳态干扰磁场和电子设备零点常偏这两个不具有独立可观性的偏置量的合成量。

作者提出一个扩维E KF 滤波器,同时估计四元数、角速率常漂和地磁观测常偏置量。

本文进一步引入更多系统误差的影响,并基于一个完整的TAM 观测模型,利用扩维UKF 算法和双重UKF 算法给出几个更加完善的地磁姿态确定算法。

1　地磁矢量观测模型和标量校正观测模型一个完整的地磁观测模型如下所示[9]:b M ,k =(I 3×3+D k )-1(O T A BR ,k b R ,k +b Biask +εk )　(1)其中,I 3×3表示单位矩阵;b M ,k 表示TAM 观测量;O 表示正交安装误差矩阵;A BR ,k 表示星体系(B 系)相对参考系(R 系)的姿态矩阵;b R ,k 为依据地磁模型得到的参考矢量;b Biask 表示常值偏置量;εk 被认为是零均值高斯地磁观测噪声,且方差为΢k =E {εk εk T}。

TAM 标量校正观测模型为[9]:z k ≡‖b M ,k ‖2-‖b R ,k ‖2=-b M ,k T(2D k +D 2k )b M ,k +　2b TM ,k (I 3×3+D k )b Biask-‖b kBias‖2+v effk (2)其中,D k 矩阵的对角线元素表示刻度因子调试误差,其它非对角线元素表示非正交性因素;v effk 表示有效噪声,且v effk ≡2[(I 3×3+D k )b M ,k -b Biask ]εk -‖εk ‖2　(3)其统计模型可用高斯模型描述,且均值和方差分别为:μeff k=-tr(΢k)(4)(σeff k)2=νT΢kν+2(tr΢2k)(5)且ν≡(I3×3+D k)b M,k-b Bias k(6)由于式(2)、(3)和(6)中的b M,k是未知量,所以必须用真实测量值b M,k替代。

2　同时校正三轴磁强计的姿态确定算法这部分基于扩维估计方案和双重滤波[10]方案给出1个姿态确定扩维UKF滤波器(记为AUKF)和2个双重UKF滤波器(记为DUKF)。

2.1　系统方程状态量定义为:x k≡[q T k,βk T,ΘT k]T(7)其中,q k表示四元数;βk表示速率常漂;9维系统误差参数Θk定义为:Θk≡(b Bias k)T　[D k,11　D k,22　D k,33　D k,12　D k,13　D k,23]TT(8) 离散四元数运动学方程[4]与陀螺常漂及TAM系统误差参数扩展而成的系统方程如下所示:q k+1βk+1Θk+1=Υk+1,k(ωk-βk-ηυ,k)q kβk+ηu,kΘk(9)其中,ωk表示角速率测量值;ηυ,k可看成是零均值高斯白噪声,且方差阵为σ2υI3×3/Δt;ηu表示常漂量的一阶Markov模型随机驱动噪声,可看成是零均值高斯白噪声,且方差阵为σ2u I3×3Δt。

2.2　AUKF滤波器文献[4]提出一个经典姿态确定UKF算法,即USQUE。

为了解决四元数估计1度冗余性问题,US-QUE采用了约束估计思想,即依据四元数后验估计利用四元数运动学方程预报一个参考四元数,并利用一个非冗余的3分量姿态参数来表示并估计真实四元数与参考四元数之间的偏差四元数。

该文基于TAM系统误差得到完美校正的假设(即Θ=09×1)给出一个具体的地磁姿态确定算法。

但在现实中,这种假设往往是不成立的。

完整地磁矢量观测模型(1)显示,地磁矢量观测与四元数表示的姿态矩阵A BR,k以及TAM系统误差因素Θk都直接相关,因而利用地磁观测获得的新息同样可以对Θk进行观测更新估计。

(关于这个扩维状态系统的可观性分析请见3.3.1小节)需要特别指出的是,在地磁矢量为唯一方向矢量观测条件下,误装矩阵O是不可观的,而要使之成为可观,则必须引入其它高精度方向矢量观测[9]。

为了分析简单起见,本文假设O=I3×3。

通过直接将该参数扩展到原状态量中得到一个16维扩维状态量(7),其系统方程和观测方程分别为式(9)和式(1)。

基于此,依据USQUE不难得到同时估计参数矢量Θk的AUKF,其实际状态量(15维)定义为X aug k=[δP k T　βT k　ΘT k]T其中,δP k为偏差四元数对应的3维GRP姿态参数[4,15]。

然后,利用经典UKF理论就可以给出相应的AUKF,而四元素均值估计则由GRP姿态参数估计(转换成偏差四元数)与预报参考四元数合成得到[4]。

需要注意的是,每次观测更新结束之后必须将δP k部分重置为零矢量。

至于该滤波器的具体形式,这里就不再详细给出了。

2.3　DUKF滤波器DUKF包括一个姿态确定UKF滤波器和一个TAM系统误差参数UKF滤波器。

本文提出一个解耦DUKF滤波器和一个耦合DUKF滤波器。

所谓耦合或解耦指的是TAM系统误差校正与姿态信息(A BR,k(q k))有关或无关。

具体来说就是,解耦DUKF的姿态确定USQUE采用简单地磁观测模型,其参数校正UKF[8]则选择与姿态无关的标量校正观测模型(2)。

耦合DUKF的姿态确定UKF同样选择USQUE,但采用完整地磁观测模型(1),其参数校正UKF也采用同一完整观测模型。

为描述简洁起见,这里仅仅介绍一下耦合DUKF的参数UKF:(1)滤波器初始化(k=0)记初始参数均值和协方差估计分别为Θ0和PΘ0,令Θ0|0=Θ0,PΘ0|0=PΘ0。

(2)时间更新(k=1,…,N)参数均值和协方差时间更新为:Θk|k-1=Θk-1|k-1PΘk|k-1=PΘk-1|k-1 然后,依据Θk|k-1和PΘk|k-1,利用确定性UT采样[8]得到t k时刻的参数西格玛点集{Θi k|k-1}180,以及均值计算权值{W(m,i)k}18i=0和协方差计算权值{W(c,i)k}18i=0;(3)观测更新(k=1,…,N)预报观测量西格玛点集为855第3期黄　琳等:卫星姿态确定与三轴磁强计校正z ik =(I 3×3+D i k)-1[A BR ,k (q-k |k -1)b R ,k +b Bias ,ik ],i =0,…,18其中,D i k表示由Θi k |k -1最后6个分量组成的D k 矩阵;bBias ,ik为Θik |k -1前3个分量组成的偏置矢量;q-k |k -1依据前一时刻USQUE 给出的四元数后验估计q-k -1|k -1和陀螺角速率观测值利用离散四元数运动学方程(9)预报得到。

为了利用三轴陀螺(相对TAM )的高频采样信息(假设采样频率相差K RI G 倍),本文提出如下方法:q-j +1|j =Υj +1,j ( ωk -1+j -βk -1+j |k -1)q-j |j -1, i =0,…,18;j =0,…,K RI G -1当j =0时,q-j |j -1=q-k -1|k -1;当j =K R IG }-1时,令q-k |k -1=q-j +1|j ;βk -1+j |k -1为USQUE 给出的陀螺常漂估计βk -1|k -1的时间更新,且βk -1+j |k -1=βk -1|k -1; ωk -1+j 表示地磁采样间隔[t k -1,t k ]内第j 个陀螺测量值。

预报观测量均值为zk=∑18i =0W(m ,i )k zik 根据地磁粗糙测量z k = b M ,k ,分别得到新息及其协方差υΘk=z k -zkP υυ,Θk=∑24i =0W (c ,i )k(z ik -zk)(z ik -z k)T+R k其中,R k 表示观测噪声(I 3×3+D k )-1εk 的协方差阵,且R k =(I 3×3+Dk )-1΢k (I 3×3+Dk )-T其中,Dk 为先验均值估计Θk |k -1后6个分量组成的D k 矩阵;(·)-T表示逆矩阵转置操作。

相关协方差阵和增益矩阵分别为:Pxz ,Θk =∑18i =0W(c ,i )k(Θik |k -1-Θk |k -1)(z ik-zk)TK Θk =P xz ,Θk(P υυ,Θk)-1 最后,参数估计观测更新为Θ k |k =Θ k |k -1+K Θk υΘkP Θk |k =P Θk |k -1-K Θk P υυ,Θk(K Θk )T 存储Θk |k 和P Θk |k ,以备USQUE 和参数UKF 算法下次递推估计之用。