三轴磁强计观测的姿态确定方法

三轴磁强计观测的姿态确定方法

0在线校正tam系统误差三轴磁强计观测（三轴磁强计，简单记录为tam）观测的近地卫星姿态确定方法已经开发。但是,大多数方法都采用简单的地磁观测模型,仅仅考虑了随机观测噪声的影响,而没有完整考虑其它诸多系统误差因素的影响。比较典型的系统误差影响因素有TAM本身三个磁通线圈轴的非正交性、刻度因子调试误差、零点常偏,以及剩磁场干扰和电磁设备工作时产生的附加干扰磁场。一些实际案例也证实对TAM进行在线校正以确保获得一定精度的姿态确定结果来说是非常必要的。Crassidis等人提出三个适合在线校正TAM系统误差(3个非正交因素、3个刻度因子误差和3个常偏置量)的参数滤波算法,即序贯中心化方法、普通EKF和简单UKF校正算法。事实上,后两者的实质区别在于新息协方差阵的计算有所不同,但UKF滤波器通常能够取得相对其它两个滤波器来说更高一些的估计精度和更快的收敛速度,而且对有色观测噪声也显示出良好的鲁棒性。需要指出的是,这三个校正算法仅仅利用了地磁观测量的大小,而与其指向(或卫星姿态)无关,属于所谓的标量校正算法。目前,利用非线性滤波理论解决姿态确定与TAM在线校正这一组合估计问题还是一个新鲜课题。本文作者曾经讨论了地磁姿态确定与地磁观测常偏置量校正两者组合的估计问题,该偏置量实际上是稳态干扰磁场和电子设备零点常偏这两个不具有独立可观性的偏置量的合成量。作者提出一个扩维EKF滤波器,同时估计四元数、角速率常漂和地磁观测常偏置量。本文进一步引入更多系统误差的影响,并基于一个完整的TAM观测模型,利用扩维UKF算法和双重UKF算法给出几个更加完善的地磁姿态确定算法。1k33+dk的校正一个完整的地磁观测模型如下所示:bΜ,k=(Ι3×3+Dk)-1(ΟΤABR,kbR,k+bBiask+εk)(1)其中,I3×3表示单位矩阵;bM,k表示TAM观测量;O表示正交安装误差矩阵;ABR,k表示星体系(B系)相对参考系(R系)的姿态矩阵;bR,k为依据地磁模型得到的参考矢量;bBiask表示常值偏置量;εk被认为是零均值高斯地磁观测噪声,且方差为Σk=E{εkεkT}。TAM标量校正观测模型为:zk≡∥bΜ,k∥2-∥bR,k∥2=-bΜ,kΤ(2Dk+D2k)bΜ,k+2bΤΜ,k(Ι3×3+Dk)bBiask-∥bkBias∥2+veffk(2)其中,Dk矩阵的对角线元素表示刻度因子调试误差,其它非对角线元素表示非正交性因素;veffk表示有效噪声,且veffk≡2[(Ι3×3+Dk)bΜ,k-bBiask]εk-∥εk∥2(3)其统计模型可用高斯模型描述,且均值和方差分别为:μeffk=-tr(Σk)(4)(σeffk)2=νΤΣkν+2(trΣ2k)(5)且ν≡(Ι3×3+Dk)bΜ,k-bBiask(6)由于式(2)、(3)和(6)中的bM,k是未知量,所以必须用真实测量值˜bM,k替代。2同时，修正三轴磁强计的状态这部分基于扩维估计方案和双重滤波方案给出1个姿态确定扩维UKF滤波器(记为AUKF)和2个双重UKF滤波器(记为DUKF)。2.1dk模型和tam系统误差参数状态量定义为:xk≡[qΤk,βkΤ,ΘΤk]Τ(7)其中,qk表示四元数;βk表示速率常漂;9维系统误差参数Θk定义为:Θk≡[(bBiask)Τ[Dk,11Dk,22Dk,33Dk,12Dk,13Dk,23]Τ]Τ(8)离散四元数运动学方程与陀螺常漂及TAM系统误差参数扩展而成的系统方程如下所示:[qk+1βk+1Θk+1]=[Φk+1,k(˜ωk-βk-ηυ,k)qkβk+ηu,kΘk](9)其中,˜ωk表示角速率测量值;ηυ,k可看成是零均值高斯白噪声,且方差阵为σ2υI3×3/Δt;ηu表示常漂量的一阶Markov模型随机驱动噪声,可看成是零均值高斯白噪声,且方差阵为σ2uI3×3Δt。2.2基于地磁矢量观测的k文献提出一个经典姿态确定UKF算法,即USQUE。为了解决四元数估计1度冗余性问题,USQUE采用了约束估计思想,即依据四元数后验估计利用四元数运动学方程预报一个参考四元数,并利用一个非冗余的3分量姿态参数来表示并估计真实四元数与参考四元数之间的偏差四元数。该文基于TAM系统误差得到完美校正的假设(即Θ=09×1)给出一个具体的地磁姿态确定算法。但在现实中,这种假设往往是不成立的。完整地磁矢量观测模型(1)显示,地磁矢量观测与四元数表示的姿态矩阵ABR,k以及TAM系统误差因素Θk都直接相关,因而利用地磁观测获得的新息同样可以对Θk进行观测更新估计。(关于这个扩维状态系统的可观性分析请见3.3.1小节)需要特别指出的是,在地磁矢量为唯一方向矢量观测条件下,误装矩阵O是不可观的,而要使之成为可观,则必须引入其它高精度方向矢量观测。为了分析简单起见,本文假设O=I3×3。通过直接将该参数扩展到原状态量中得到一个16维扩维状态量(7),其系统方程和观测方程分别为式(9)和式(1)。基于此,依据USQUE不难得到同时估计参数矢量Θk的AUKF,其实际状态量(15维)定义为Xaugk=[δΡkΤβΤkΘΤk]Τ其中,δPk为偏差四元数对应的3维GRP姿态参数。然后,利用经典UKF理论就可以给出相应的AUKF,而四元素均值估计则由GRP姿态参数估计(转换成偏差四元数)与预报参考四元数合成得到。需要注意的是,每次观测更新结束之后必须将δPk部分重置为零矢量。至于该滤波器的具体形式,这里就不再详细给出了。2.3耦合或解耦的姿态确定模型DUKF包括一个姿态确定UKF滤波器和一个TAM系统误差参数UKF滤波器。本文提出一个解耦DUKF滤波器和一个耦合DUKF滤波器。所谓耦合或解耦指的是TAM系统误差校正与姿态信息(ABR,k(qk))有关或无关。具体来说就是,解耦DUKF的姿态确定USQUE采用简单地磁观测模型,其参数校正UKF则选择与姿态无关的标量校正观测模型(2)。耦合DUKF的姿态确定UKF同样选择USQUE,但采用完整地磁观测模型(1),其参数校正UKF也采用同一完整观测模型。为描述简洁起见,这里仅仅介绍一下耦合DUKF的参数UKF:(1)滤波器初始化k.0记初始参数均值和协方差估计分别为Θ0和PΘ0,令ˆΘ0|0=Θ0‚ΡΘ0|0=ΡΘ0。(2)确定参数下的西格玛点集参数均值和协方差时间更新为:ˆΘk|k-1=ˆΘk-1|k-1ΡΘk|k-1=ΡΘk-1|k-1然后,依据ˆΘk|k-1和PΘk|k-1,利用确定性UT采样得到tk时刻的参数西格玛点集{Θik|k-1}180,以及均值计算权值{W(m,i)k}18i=0和协方差计算权值{W(c,i)k}18i=0;(3)方法及估计结果预报观测量西格玛点集为zik=(Ι3×3+Dik)-1[ABR,k(ˆq-k|k-1)bR,k+bBias,ik],i=0,⋯,18其中,Dik表示由Θik|k-1最后6个分量组成的Dk矩阵;bBias,ik为Θik|k-1前3个分量组成的偏置矢量;ˆq依据前一时刻USQUE给出的四元数后验估计ˆq¯k-1|k-1和陀螺角速率观测值利用离散四元数运动学方程(9)预报得到。为了利用三轴陀螺(相对TAM)的高频采样信息(假设采样频率相差KRIG倍),本文提出如下方法:ˆq¯j+1|j=Φj+1,j(˜ωk-1+j-ˆβk-1+j|k-1)ˆq当j=0时,ˆq;当j=KRIG}-1时,令ˆq;ˆβk-1+j|k-1为USQUE给出的陀螺常漂估计ˆβk-1|k-1的时间更新,且ˆβk-1+j|k-1=ˆβk-1|k-1;˜ωk-1+j表示地磁采样间隔[tk-1,tk]内第j个陀螺测量值。预报观测量均值为ˆzk=18∑i=0W(m,i)kzik根据地磁粗糙测量zk=˜bΜ,k,分别得到新息及其协方差υΘk=zk-ˆzkΡυυ,Θk=24∑i=0W(c,i)k(zik-ˆzk)(zik-ˆzk)Τ+Rk其中,Rk表示观测噪声(I3×3+Dk)-1εk的协方差阵,且Rk=(Ι3×3+ˆDk)-1Σk(Ι3×3+ˆDk)-Τ其中,ˆDk为先验均值估计ˆΘk|k-1后6个分量组成的Dk矩阵;(·)-T表示逆矩阵转置操作。相关协方差阵和增益矩阵分别为:Ρxz,Θk=18∑i=0W(c,i)k(Θik|k-1-ˆΘk|k-1)(zik-ˆzk)ΤΚΘk=Ρxz,Θk(Ρυυ,Θk)-1最后,参数估计观测更新为ˆΘk|k=ˆΘk|k-1+ΚΘkυΘkΡΘk|k=ΡΘk|k-1-ΚΘkΡυυ,Θk(ΚΘk)Τ存储ˆΘk|k和PΘk|k,以备USQUE和参数UKF算法下次递推估计之用。3计算机模拟和结果分析3.1真实地磁场仿真本文考虑一个典型的近地刚体小卫星。其近圆轨道(高度为823km、倾角为82deg)数据由商业软件STK╋产生。卫星处于无控状态,但受重力梯度力矩、剩磁力矩(剩磁矩为0.3A·m2)等内/外界干扰力矩的影响,其初始旋转角速率为2.0deg/s。真实地磁观测量由真实地磁场矢量加上各种干扰因素的影响合成得到。本文采用一个10×10阶的IGRF模型仿真出真实地磁场。各种具体干扰因素综合形成的系统误差影响假设为:bBiastrue=[5000nΤ3000nΤ6000nΤ]ΤDtrue=[0.050.050.050.050.100.050.050.050.05]另外,地磁矢量观测还受到随机观测噪声εk的影响,具体有TAM敏感器噪声、地磁场模型化误差以及定轨误差等。实际观测噪声εk是一个有色噪声,本文利用一阶Markov模型来产生这一观测噪声,其时间常数对应于18deg轨道弧长(大概300s左右),有色噪声大小与标准差为60nT的高斯白噪声相同,采样周期为10s。三轴陀螺真实角速率测量值由卫星真实旋转角速率以及常漂参数和随机高斯噪声的影响合成得到。真实常漂为0.1deg/h,而σ2υ=0.1(μrad)2/s,σu2=1×10-7(μrad)2/s3,采样周期为1s。3.2启动滤波和参数调整(1)初始协方差阵的选择仿真中选择的初始四元数估计与真实四元数之间的偏差四元数的欧拉角大概为50deg左右。令δP0=03×1,其初始协方差阵则选择一个相当大的矩阵(尽管不具有物理意义,但能够加快姿态估计收敛的速度)。陀螺常漂初始估计为0.2deg/h,其初始协方差估计选为Ρ0β=diag([0.2deg/h0.2deg/h0.2deg/h]2)(2)差阵选设置令初始参数均值估计Θ0=09×1;其协方差阵选为Ρ0Θ=diag([1e-5nΤ1e-5nΤ1e-5nΤ0.10.10.10.10.10.1]2)(3)参考地磁变量的生成与滤波器观测噪声的选择利用已知轨道信息和一个8×8阶IGRF模型产生所需的参考地磁矢量。滤波器采用的地磁观测噪声标准差经过调试选为5×60nT(σ)。3.3模拟结果与分析3.3.1不同线性化系统及的异质性这部分考察了基于地磁矢量观测估计四元数、陀螺常漂参数以及TAM校正参数这一扩维状态的可观性。连续四元数运动学方程与陀螺常漂及TAM系统误差参数扩展而成的状态方程和观测方程(1)构成的系统是明显的非线性系统,而目前分析非线性系统可观性的技术和方法还很不成熟,一个通用的处理方法是将连续非线性系统离散线性化,然后利用线性系统可观性分析手段对非线性系统的局部可观性能进行分析。本文考虑的扩维系统经过离散线性化处理后的形式如下:Δxk+1=Ψ(tk+1,tk)Δxk+ζk,k=0,⋯,Ν-1(10)zk+1=ΗkΔxk+ξk,k=1,⋯‚Ν(11)式中,Δxk≡[ΔqkΤ,ΔβkΤ,ΔΘkΤ]T,表示真实状态与估计状态之间的误差,其中,Δq表示四元数估计误差的矢量部分;ζk和ξk分别表示离散后的系统噪声和观测噪声;Ψ(tk+1,tk)表示状态转移矩阵,且Ψ(tk+1,tk)≈I15×15+FkΔt,其中,Δt表示离散间隔,I15×15表示单位矩阵,而线性化系统矩阵Fk不难从文献得知:Fk=[-[ω˜k×]3×3-12Ι3×303×9012×3012×3012×9]式(11)中,Hk表示观测矩阵,根据文献,并经过一定合理简化后得到:Ηk=[2[B^ref×]3×303×3Ι3×3-Σ]上式中,B^ref=ABR,k(q^¯k)bR,k,而Σ=[Br,100Br,2Br,30Br,20Br,10Br,3Br,30Br,1Br,2]其中,Br=[Br,1Br,2Br,3]T表示参考地磁矢量观测值。上述局部线性化离散系统(10)和(11)的可观性可根据一段观测时间区间上的Gramian矩阵的奇异性与否判断得知,ϑ=∑k=1ΝΨΤ(tk,t0)ΗkΤΗkΨ(tk,t0)若上述矩阵是非奇异矩阵,则该系统是可观的。令Δt=10s,N=3600,通过研究发现,各种已考虑到情况下的Gramian矩阵都是非奇异的,证实该扩维状态系统确实是可观的。另外,通过对该矩阵的特征值分解,发现其最小特征值与最大特征值的比率约为3×10-10;若忽略掉9维TAM待校正参数相应的行和列,发现剩余6维子系统矩阵的最小特征值与最大特征值之间的比率增大为5×10-6,这表明TAM待校正参数的可观性相对其它6维姿态参数来说要弱一些。3.3.2不同滤波器下的姿态确定误差这部分考察了经典USQUE(未校正TAM系统误差)、AUKF和解耦/耦合DUKF共4个姿态确定滤波器在地磁矢量观测存在系统误差和有色观测噪声情况下的姿态确定性能。本文利用文献提供的姿态确定性能评估方法,即,利用真实四元数与估计四元数之间的偏差四元数的旋转欧拉角(单位为deg)评估四元数估计好坏,δα=2arccos(δq4)其中,δq4表示标量部分;利用陀螺常漂估计β^Est与真实值βTrue之间的误差的模(单位deg/h)评估其估计好坏|Δβ|=|β^Est-βΤrue|图1和图2给出4个滤波器的一组典型姿态确定误差结果。显然,未校正TAM观测系统误差的USQUE的姿态确定误差远大于AUKF、解耦DUKF以及耦合DUKF的估计误差。现实中,由于TAM观测无可避免地存在某种程度的系统误差或出现某些极端情况,因而校正TAM观测量对于保障姿态确定算法给出一定精度的姿态确定结果来说意义重大。这部分还利用蒙特卡洛仿真方法比较了AUKF、解耦DUKF及耦合DUKF三者的估计性能。在20次独立仿真过程中,这三者都取得了明显收敛的结果。图3、图4分别给出了TAM待校正参数(bBias的第三个分量和D33分量)估计误差的均方根RMS(RootMeanSquare)结果。不难看出,这三个算法之所以能够取得远远好于经典USQUE算法姿态确定性能