本发明属于控制领域，涉及一种航天器的联合姿态确定方法。
背景技术
为了实现高精度地球重力场信息测量，重力梯度测量卫星携带的有效载荷重力梯度仪具有极高的测量精度。为了避免对高精度有效载荷的工作性能造成破坏，通常使用的三浮陀螺由于具有活动部件而不被容许使用。因此，不能采用陀螺+星敏感器实现姿态确定。考虑到重力梯度仪除了提供高精度的线加速度测量信息以外，还可以提供角加速度测量信息，因此，提出重力梯度仪及星敏感器的姿态确定算法。
技术实现要素：
本发明解决的技术问题是：在轨卫星均采用陀螺+星敏感器的方式或无陀螺模式实现姿态确定，考虑到重力梯度测量卫星不能携带具有活动部件的三浮陀螺，因此提出了在平台不具备角速度测量单元条件下的一种新型姿态确定方法，即角加速度测量元件+星敏感器的联合定姿算法。
本发明的技术解决方案是：
一种重力梯度仪及星敏感器的联合定姿方法，步骤如下：
(1)建立重力梯度仪的测量方程，得到重力梯度仪的观测量；
(2)建立离散化的状态方程；
(3)根据所述离散化的状态方程，建立状态预测误差的相关矩阵的时间更新关系；
(4)建立星敏感器的测量方程，得到星敏感器的观测量；
(5)根据步骤(2)获得的离散化的状态方程以及步骤(4)获得的观测量，确定一步滤波估计误差相关矩阵；
(6)根据步骤(1)得到的重力梯度仪观测量以及步骤(3)建立的状态预测误差的相关矩阵的时间更新关系，获得姿态、姿态角速度以及单轴重力梯度仪漂移的一步预测估计值；
(7)根据步骤(4)得到的星敏感器的观测量获得星敏感器的测量新息；
(8)根据步骤(7)得到的测量新息、步骤(6)得到的姿态、姿态角速度以及单轴重力梯度仪漂移的一步预测估计值，得到姿态四元素、姿态角速度及单轴重力梯度仪漂移的滤波估计值，从而完成重力梯度仪及星敏感器的联合定姿。
所述步骤(1)建立重力梯度仪的测量方程，得到重力梯度仪的观测量，具体为：所述重力梯度仪的观测量包括角加速度测量输出向量。
所述步骤(2)建立离散化的状态方程，具体为：
在重力梯度仪采样周期[tk-1,tk]的时间区间δt＝tk-tk-1内，建立离散化的状态方程δxk＝φk,k-1δxk-1+wk-1；
其中，φk,k-1≈i9×9+f(tk-1)δt，
所述步骤(3)建立状态预测误差的相关矩阵的时间更新关系，具体为：
其中，pk/k-1表示预测状态误差的相关矩阵，pk-1表示t＝tk-1时刻的一步滤波估计误差相关矩阵。
所述步骤(4)建立星敏感器的测量方程，得到星敏感器的观测量，具体为：
对于星敏感器a和星敏感器b，建立星敏感器测量方程为：
其中，zsa为星敏感器a光轴的方向矢量在卫星本体坐标系中的测量值，zsab为星敏感器a光轴的方向矢量在卫星本体坐标系中坐标的标称值，为长矢量；
一步滤波估计误差相关矩阵pk，具体为：
所述步骤(6)姿态、姿态角速度以及单轴重力梯度仪漂移的一步预测估计值通过如下方式获得：
其中，
测量新息为δz＝zsa+zsb。
姿态四元素、姿态角速度及单轴重力梯度仪漂移的滤波估计值具体通过如下方式获取：
其中
本发明与现有技术相比的优点在于：
(1)本发明提出的重力梯度仪及星敏感器联合定姿算法，充分利用有效载荷提供的角加速度信息，建立包括角加速度在内的三阶动力学方程，实现包括姿态角、姿态角速度及姿态角加速度的预测，这与以往仅考虑姿态角速度及姿态角的二阶动力学方程进行姿态预测的方法相比，具有很大的创新性。
(2)本发明提出的重力梯度仪及星敏感器联合定姿算法，可推广应用于没有角速度测量敏感器，仅有角加速度测量敏感器及姿态角测量敏感器卫星平台的姿态确定。
附图说明
图1为本发明方法流程图；
图2为三轴姿态曲线示意图。
具体实施方式
目前，在轨卫星均采用陀螺+星敏感器的方式或无陀螺模式实现姿态确定，考虑到重力梯度测量卫星不能携带具有活动部件的三浮陀螺，因此提出了在平台不具备角速度测量单元条件下的一种新型姿态确定方法，即一种重力梯度仪及星敏感器的联合定姿方法，如图1所示，步骤如下：
(1)建立重力梯度仪的测量方程，得到重力梯度仪的观测量；具体为：所述重力梯度仪的观测量包括角加速度测量输出向量。
角加速度测量模型(由重力梯度仪差模测量模型提供)
y(t)＝a(t)+b(t)+η1(t)
式中y(t)是角加速度测量输出向量；a(t)为重力梯度卫星实际的角加速度；b(t)是重力梯度仪测量模型中的常值漂移；η1(t)、η2(t)均近似为三维高斯白噪声向量，即：
e[η1(t)]＝0，
e[η2(t)]＝0，
(2)建立离散化的状态方程；
具体为：
首先设t时刻的四元数、姿态角速度和重力梯度仪差模测量模型中的常漂估计值为
由重力梯度仪测量值确定的角速度增量为：
其中角加速度估计值为
加速度计测量的角加速度(重力梯度仪差模)经过积分得到角速度，假设在tk-1时刻确定的姿态角速度为
进而可以求出从tk-1时刻到tk时刻姿态角机动四元数参数为：
由四元数乘法可以计算得出由tk-1时刻到tk时刻的姿态四元素状态预测为：
其中
当时间区间δt较短，姿态变化较小时，可以近似为：
相应的姿态预测由如下公式确定：
根据重力梯度仪测量值确定的角速度增量公式，姿态角速度状态预测为：
常漂的预测估计由常漂表达式可以得出：
定义如下的状态偏差向量
其中
上式中符号
其中
w(t)为观测噪声，满足：
e[w(t)]＝06×1
e[w(t)wt(t')]＝q(t)δ(t-t')
则在重力梯度仪采样周期[tk-1,tk]的时间区间δt＝tk-tk-1内，离散化的状态方程为：
δxk＝φk,k-1δxk-1+wk-1
上式中
qk≈q(t)·δt
其中，φk,k-1≈i9×9+f(tk-1)δt，
(3)根据所述离散化的状态方程，建立状态预测误差的相关矩阵的时间更新关系；
具体为：
其中，pk/k-1表示预测状态误差的相关矩阵，pk-1表示t＝tk-1时刻的一步滤波估计误差相关矩阵。
(4)建立星敏感器的测量方程，得到星敏感器的观测量；
具体为：
对于星敏感器a和星敏感器b，建立星敏感器测量方程为：
其中，zsa为星敏感器a光轴的方向矢量在卫星本体坐标系中的测量值，zsab为星敏感器a光轴的方向矢量在卫星本体坐标系中坐标的标称值，为长矢量；
(5)根据步骤(2)获得的离散化的状态方程以及步骤(4)获得的观测量，确定一步滤波估计误差相关矩阵；具体为：
根据步骤(1)得到的重力梯度仪观测量以及步骤(3)建立的状态预测误差的相关矩阵的时间更新关系，获得姿态、姿态角速度以及单轴重力梯度仪漂移的一步预测估计值；
其中，
(7)根据步骤(4)得到的星敏感器的观测量获得星敏感器的测量新息δz＝zsa+zsb；
(8)根据步骤(7)得到的测量新息、步骤(6)得到的姿态、姿态角速度以及单轴重力梯度仪漂移的一步预测估计值，得到姿态四元素、姿态角速度及单轴重力梯度仪漂移的滤波估计值，从而完成重力梯度仪及星敏感器的联合定姿。
姿态四元素、姿态角速度及单轴重力梯度仪漂移的滤波估计值具体通过如下方式获取：
其中
实施例：
设定动力学仿真步长为25ms，控制器采样及控制周期为0.1s，取星敏感器光轴/横轴测量精度(3σ)分别为3″/35″，重力梯度仪测量噪声在测量带宽内为
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。