一种用于卫星姿态确定的交互式滤波方法
【技术领域】：
[0001] 本发明涉及一种用于卫星姿态确定的交互式滤波方法，属于利用最优K矩阵方法 和非线性滤波技术进行卫星姿态估计的技术领域。
【背景技术】：
[0002] 卫星姿态估计技术是航天技术的关键技术之一，由陀螺与星敏感器组成的卫星姿 态估计系统由于测姿精度高、可靠性好以及自主性强等优点得到了广泛的应用。K矩阵姿 态估计方法是基于Wahba问题的四元数确定方法，具有计算效率高的优点。在QUEST算法 的基础上，Itzhack Y. Bar-Itzhack提出了利用陀螺量测更新K矩阵的REQUEST算法，进一 步提高了 QUEST算法的估计精度，但是REQUEST算法具有增益因子不可调节的缺点，无法 根据实际的量测误差，自动调节增益因子，使得算法的适用性存在一定问题。在此基础上， D. Choukroun等于2004年提出了 optimal-REQUEST算法，在经典Kalman滤波框架的基础 下，利用最优增益因子提高了姿态估计精度。但是optimal-REQUEST算法存在不能估计陀 螺常值漂移的缺点，而且最优增益因子无法补偿陀螺常值漂移带来的影响，使得算法存在 发散现象。
[0003] 为克服optimal-REQUEST算法无法估计陀螺常值漂移的缺陷，采用CKF算法与 optimal-REQUEST算法相结合的方式实现陀螺常值漂移的估计，提高算法的适用性。由于陀 螺常值漂移估计为非线性状态空间模型，常用的方法有EKF、UKF和CKF算法，但是EKF算法 存在理论上的局限性：模型线性化引入了截断误差，导致滤波精度下降，同时要计算雅克比 矩阵，计算复杂。UKF与CKF算法同是基于最优高斯滤波框架，相比于UKF，CKF采用三阶球面 相径容积规则来近似非线性函数的均值和方差，可以保证在理论上以三阶多项式逼近任何 非线性高斯状态的后验均值和方差，具有实现简单，收敛性好等优点。将optimal-REQUEST 算法与CKF算法相结合的滤波方式，可以发挥各自的优势，optimal-REQUEST算法具有线性 特征，估计最优四元数，可以避免四元数规范化带来的误差，而采用CKF算法估计陀螺常值 漂移，也可以解决量测方程非线性化带来的估计精度差的缺点，使得算法的估计精度得到 提高，提高了算法的适用性。
【发明内容】
：
[0004] 本发明的目的是为了提高姿态估计精度和增强系统的适用性，提出了一种用于 卫星姿态确定的optimal-REQUEST与CKF交互式滤波方法。发明采用optimal-REQUEST 算法估计卫星姿态四元数，采用CKF算法估计陀螺常值漂移，并将optimal-REQUEST算 法估计的最优四元数，作为CKF算法的输入，同时CKF算法估计的陀螺常值漂移作为 optimal-REQUEST算法的更新补偿，实现两种算法的交互，达到系统的最优性能。
[0005] 本发明的技术方案具体如下：
[0006] 本发明的用于一种用于卫星姿态确定的交互式滤波方法，包括以下几个步骤：
[0007] 步骤1 :读取传感器量测数据，所述传感器测量数据包括陀螺仪数据和星敏感器 数据；
[0008] 步骤2 :建立卫星姿态估计系统的状态空间模型；
[0009] 所述的卫星姿态估计系统状态空间模型包括：
[0010] (1)陀螺常值漂移估计状态空间模型
[0015] 式中：β k表示k时刻陀螺常值漂移；β k+1为k+Ι时刻陀螺常值漂移；ij4表示k时 刻最优估计四元数；At表示采样时间；Iik为零均值高斯白噪声，方差为
为四 元数更新矩阵；Ok表示k时刻载体转速，[ω kX]表示向量叉乘矩阵；rk+1为k+Ι时刻向量 观测器参考向量；bk+1为k+Ι时刻向量观测器量测向量；氣为陀螺仪量测值，量测噪声方差 为Sbk+1为零均值高斯白噪声，方差为
'表示姿态旋转矩阵。以上 是单向量观测量测模型，对于一般星敏感器，在同一时间往往能够观测多颗有效星，其量测 模型可以扩充为：
[0017] 式中，Zk+lil表示第i个有效矢量观测；Sbk+lil为相应的随机噪声；N(k+1)表示k+1 时刻的有效向量个数；X k表示组合之后的状态参量，此处表示陀螺常值漂移β k;vk表示组 合之后的零均值高斯白噪声；h( ·)表示非线性量测函数；
[0018] (2)基于optimal-REQUEST方法的姿态确定状态空间模型
[0026] 式中：\表示状态K矩阵；Y k+1为k+1时刻量测矩阵；V k+1表示k+1时刻量测噪声； bk+i表示k+Ι时刻第i个有效观测向量；r k+lii表示k+Ι时刻第i个有效参考向量；tr ( ·) 表示取矩阵的迹；ak+lil表示k+Ι时刻第i个加权系数;WkS系统过程噪声矩阵；ε为陀螺 量测随机噪声向量；
[0027] 步骤3 :针对上述状态空间模型，在已知时刻陀螺常值漂移的状态估计和估计方 差的情况下，利用标准容积卡尔曼滤波进行陀螺常值漂移估计的时间更新和量测更新，得 到陀螺常值漂移的一步状态预测方差、量测预测方差以及互协方差；
[0028] 步骤3.1时间更新；
[0029] 已知k时刻的状态估计为爲&以及估计协方差Pklk，求取容积点为：
[0032] 式中，η为状态变量的维数；ξ i为第i个容积点；e i表示η维单位列向量在第i个 位置为I ;Pk|k表示k时刻估计状态协方差矩阵；表示k时刻状态估计；a lik|k表示传递 容积点；
[0033] 由于陀螺常值漂移状态方程是线性方程，故容积点经过状态传递之后仍为原值； 一步状态预测估计与估计方差：
[0036] 式中，表示一步预测状态；Pk+1|k表示一步预测状态协方差矩阵；Qk表示k时刻 过程噪声方差阵；
[0037] 步骤3. 2量测更新；
[0038] 计算更新容积点：
[0040] 容积点经过量测非线性函数传递为：
[0041] T i,k+i|k= h ( λ i k+1|k)
[0042] 式中，h( ·)的定义如前向量观测器非线性量测模型；λ i k+1|k为计算更新容积点；
[0043] 量测一步预测、预测方差和互协方差为：
[0047] 式中，矣+:ψ表示k+Ι时刻量测一步预测；PZZik+1|k表示k+Ι时刻量测一步预测协方差 矩阵；PXZlk+i|k表示k+Ι时刻量测与状态一步预测互协方差矩阵；Rk+1表示k+Ι时刻量测噪声 方差阵；
[0048] 陀螺常值漂移估计状态空间滤波增益、k+Ι时刻的状态估计及估计方差如下式所 示：
[0050] 式中，％+1表示k+Ι时刻容积卡尔曼滤波增益矩阵；處+1|M表示k+Ι时刻陀螺常值 漂移最优估计；尾+1表示k+Ι时刻向量观测器量测向量；Pk+1|k+1表示k+Ι时刻估计状态误差 协方差矩阵；
[0051] 步骤4 :利用步骤3估计的陀螺常值漂移补偿陀螺量测，然后利用补偿后的陀螺量 测构造四元数更新矩阵，在optimal-REQUEST算法更新过程中，修正由于陀螺常值漂移引 起的O k不准确问题，在利用optimal-REQUEST算法实现最优K矩阵的求取：
[0052] 步骤 4. I optimal-REQUEST 时间更新；
[0053] 已知k时刻的状态K矩阵估计为矣以及估计协方差％，得到一步状态预测估计 和估计方差：
[0056] 式中，i>k表示k时刻四元数更新矩阵；表示k时刻状态K矩阵一步预测； 表示k时刻一步预测状态协方差矩阵；必表示k时刻过程噪声方差阵；
[0057] 步骤 4. 2 optimal-REQ