卫星三轴姿态和姿态角速率的确定

卫星三轴姿态和姿态角速率的确定

1基本粒子和金相结合的定姿算法近年来，小型卫星技术是一个重要的指导概念，旨在降低成本。因此,仅仅使用磁强计完成三轴姿态、姿态角速率的确定无疑有着很好的应用前景。在以往的卫星姿态确定系统中,磁强计往往做为辅助手段与其它部件配合使用,常见的处理方法多为单点算法确定姿态(例如双矢量定姿),而角速率的获取需借助于陀螺仪。本文的定姿算法仅利用三轴磁强计的测量数据,通过卡尔曼滤波能同时给出三轴的姿态角和姿态角速率,从而解决了没有陀螺的情况下角速率获取问题。算法中的一些重要公式在以往的基础上做了改进,特别是作为算法核心的测量矩阵有着独到之处。此系统适用于有一定倾角的近地轨道,且卫星的定姿精度要求不高。它可作为卫星姿态控制系统中闭合回路的敏感器定姿部分。2模型的数学建模2.1卫星数学模型本文所考虑的卫星模型是重力梯度稳定,三轴主动磁控,俯仰轴负向安装偏置动量轮的卫星。卫星轨道为具有一定倾角的近地轨道,从而保证有丰富的地磁资源和三轴姿态的可观测性。下面给出的动力学和运动学方程作为滤波器的状态转移方程。对于其它的卫星模型,相应建立其动力学模型后,此算法同样适用。Ι˙ω+∼ω(Ιω+hw)=∑Τ(1)[˙q1˙q2˙q3˙q4]=12[q4-q3q2q1q3q4-q1q2-q2q1q4q3-q1-q2-q3q4][ωxωy+ω0ωz0](2)Iω˙+ω∼(Iω+hw)=∑T(1)⎡⎣⎢⎢⎢⎢q˙1q˙2q˙3q˙4⎤⎦⎥⎥⎥⎥=12⎡⎣⎢⎢⎢⎢q4q3−q2−q1−q3q4q1−q2q2−q1q4−q3q1q2q3q4⎤⎦⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢ωxωy+ω0ωz0⎤⎦⎥⎥⎥⎥(2)其中I为星体的惯量矩阵,ω为星体的绝对角速度向量,(~)表示矢量外积的矩阵形式,hw=[0-h0]T,为飞轮角动量,∑T为作用在星体上的力矩之和。q为从轨道坐标系到本体坐标系的姿态四元素,ω0为轨道角速率。2.2磁强测量值bBb=ABorb(3)Bb=ABorb(3)Bb为地磁场在星体坐标下的表示,理想情况的磁强计测量值等于Bb。Borb为地磁场在轨道坐标下的表示。A为轨道坐标系到本体坐标系的转换矩阵。2.3小卫星轨道系统的运行角速度分量和规则ΙΔ˙ω+Δ∼ω(Ιωorb+hw)+∼ωorb(ΙΔω)=ΔΤ(4)˙q=12Δω(5)Bb=[Ι-2∼q]Borb(6)其中Δω为星体相对于轨道坐标系的摄动角速度矢量,ωorb为小卫星轨道运行角速度,ΔT=∑T为作用在星体上的力矩之和,q=[q1q2q3]T。考虑星体相对于轨道坐标系的摄动为小转角(欧拉转角)情况,在线性化过程中认为q4≈1。2.4状态空间形式的离散化选取6维的状态矢量和3维的测量矢量分别为:{X=[ΔωΤqΤ]ΤY=∼BbBorb∥Bb∥∥Borb∥(7)把式(4)、(5)、(6)表示为状态空间形式,并采用MATLAB函数进行离散化:{Xk+1=AXk+BΤkYk+1=Ck+1Xk+1(8)其中Ck+1为测量矩阵,它为时变的,其公式在稍后介绍。3加密度估计卡尔曼滤波器的滤波公式和误差方差计算公式基本采用经典的算法,但在第一次估值时我们采用另外一组公式,其中Κk+1=Ρk+1CΤk+1R-1Ρ-1k+1=(Ρ´k+1)-1+CΤk+1RCk+1(9)其它公式不变。当对初始状态缺乏必要的知识时,可选取Ρ0=1εΙ(ε为很小的数)。此形式在经典公式中不便直接计算增益矩阵。构造新息序列的原则是要体现出预测值与实测值的差别,以便对状态向量的估计进行修正。由于没有任何的姿态信息表现在二者“大小”上,所以不用传统的二者之差而取二者的单位矢量的外积作为新息。E(tk+1)=∼Bmeas(tk+1)ˆB-(tk+1)∥Bmeas(tk+1)∥∥ˆB-(tk+1)∥(10)其中Bmeas为磁强计测量值,ˆB-为磁强计测量值的预测估值。得到新息序列后,对状态矢量的预测估值求偏导即可得测量矩阵,推导过程从略。C(tk+1)=-∂E(tk+1)∂ˆX-(tk+1)=[03×32∼Bmeas(tk+1)∼Borb(tk+1)Τ∥Bmeas(tk+1)∥∥ˆB-(tk+1)∥](11)四元素的更新方程不用传统的加法更新而采用四元素乘法进行更新。显然,这样做的物理意义更加明确。ˆqk+1=Δqk+1⋅q´k+1(12)q′k+1为预测的机动参数,Δqk+1为更新的机动参数。4仿真结果和分析数值仿真的思路是:假定小卫星的初始姿态ω0、q0,用解微分方程的算法求解动力学和运动学方程,得到卫星仿真时段的姿态理论值,以此作为衡量滤波器性能的标准。为了得到滤波器需要的输入参数Borb、Bmeas,我们做两点简化:假定轨道是圆轨道;地磁场模型使用偶极子模型。图1中Tg、Ts、Td分别为重力梯度力矩、仿真力矩和干扰力矩。经过分析,Td比Tg、Ts小一个数量级以上。下面给出1组仿真结果,实线为理论值,虚线为滤波估值。由图可见,滤波器的收敛时间约为3—5分钟。姿态角收敛精度在1°以内。此滤波算法主要受两方面因素的制约:其一,初始估值的影响。分析表明,俯仰角和角速率允许较大的初始偏差,另外其收敛速率、精度和稳定性均较其他两轴好。这主要是因为俯仰方向有保持惯性指向的偏置动量轮,使得未建模的干扰力矩对此轴的影响近乎为零。其二,待估实际值的取值,也即估值范围的影响本文研究的卫星模型决定了滚动轴有着小的姿态变化,其余两轴做了接近60°姿态估值,俯仰轴做了2°/s角速率估值,结果均收敛。受到线性化时“小转角”的限制,超过90°的大姿态估值使得线性化模型不再准确,导致滤波器发散,并且是无法克服的。因此,本算法在大转角的入轨捕获阶段不适用。5q合并q与r间折的关系Q、R和P-(t0)决定着滤波器的初始收敛速率和估值精度。稳态下,Q、R还决定着滤波器的稳定性。把Q/R或P-(t0)/R取得大些,有助于改善滤波器的初始收敛速率,但会导致稳态估计精度的下降,设计时必须在二者之间折中。合理的选择Q能够部分抵消模型不准和计算误差的影响,而对R的选择,包括着对磁力矩器干扰和星上剩磁干扰影响的考虑。Q取的过小或者R取的过大,都会使得滤波估值和测量序列逐渐分离而形成数据饱和。相反,模型误差逐渐积累,引起的破坏作用就越来越突出,滤波器在进入稳态跟踪后会出现较严重的偏离现象。此时,须反复调整Q,使滤波器在稳态时能保持稳定,并在规定误差范围内工作。对于一个N维状态模型,整个滤波和误差方差计算所需要的有效阶数为12Ν(Ν+3)。因此即使N增加一个很小的数目计算量也会显著增加。由于星上计算机系统总希望