CN102879011A

CN102879011A - 一种基于星敏感器辅助的月面惯导对准方法

Info

Publication number: CN102879011A
Application number: CN2012103522288A
Authority: CN
Inventors: 李骥; 张晓文; 王鹏基; 何英姿; 张洪华
Original assignee: Beijing Institute of Control Engineering
Current assignee: Beijing Institute of Control Engineering
Priority date: 2012-09-21
Filing date: 2012-09-21
Publication date: 2013-01-16
Anticipated expiration: 2032-09-21
Also published as: CN102879011B

Abstract

一种基于星敏感器辅助的月面惯导对准方法，针对的对象是装备有捷联惯导系统的月面探测器。对准的基本原理是利用加速度计确定重力方向建立水平面；利用星敏感器惯性姿态测量结合已知的起飞点位置确定方位角。在方法实现上分成了星敏+陀螺的惯性姿态估计和利用静基座导航速度误差的惯导平台修正两个步骤。本发明方法成功解决了月球表面惯导对准的问题，对准误差不大于0.05°。

Description

一种基于星敏感器辅助的月面惯导对准方法

技术领域

本发明涉及一种月球探测自主导航方法。

背景技术

对于月球采样返回这类任务来说，当完成月面采样之后，登月探测器需要携带月壤样品从月面起飞上升，完成与留轨飞行器的交会对接，并将样品交由留轨飞行器送返地球。因此月面起飞上升是月球采样返回任务成功的关键环节之一。从月球上起飞或者发射航天器采用的基本手段与地球火箭一样是惯性导航。惯性导航是一种递推式的导航方法，它需要精确的初值，因此在起飞前必须完成惯导系统的对准。地球上进行惯导对准的基本方法是自对准，即利用加速度计获得的重力方向建立水平基准，利用陀螺获得的地球自转角速度方向确定方位基准，必要时还可以通过光学瞄准设备进一步提高方位对准精度。在月球上，由于月球自转角速度很小，重力加速度也只有地球的1/6，因此惯导自对准的精度很低，特别是方位误差会达到几度甚至十几度。而且月球上没有类似地球火箭发射的光学瞄准设备，也就无法提高方位对准精度。这些因素导致在月球上进行惯导对准不能直接采用地球上的常规方法。

目前为止，国外成功完成月面起飞的只有前苏联的“月球”系列无人探测器和美国的“阿波罗”登月飞船。前苏联的“月球”探测器由于采用的是垂直上升的飞行轨迹，它只需要利用加速度计确定铅垂线方向，而无需确定方位角，因此不构成完整的惯导对准。而美国的“阿波罗”登月飞船则采用的是手动光学瞄准技术。该技术利用一种光学瞄准望远镜来实施。首先由宇航员手动瞄准某颗已知恒星，然后从望远镜上的刻度读出该恒星的方位信息，之后对另一颗恒星再次进行同样测量就可以解算出登月飞船的姿态，最后通过电路系统驱动IMU完成对准。对于无人的月球探测器来说，这种有人参与的方法显然不能使用，而且人工瞄准的方法精度也较低。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种星敏感器辅助的月面惯导自对准方法，针对装备有捷联惯导系统的月面探测器，利用加速度计确定重力方向建立水平面，利用星敏感器惯性姿态测量结合已知的起飞点位置确定方位角，从而解决了月球表面惯导对准的问题。

本发明的技术解决方案是：一种基于星敏感器辅助的月面惯导对准方法，步骤如下：

(1)利用陀螺测量得到探测器的角速度

扣除陀螺常值漂移得到探测器的角速度估计值

其中下标k表示对应时间为t_k，上标b表示探测器本体系；

(2)利用姿态运动学方程对探测器的惯性姿态四元数进行更新，其中q是表征探测器本体系相对惯性系的姿态四元数，且有q＝[q₁ q₂ q₃ q₄]^T，下标k|k-1表示由时间t_k-1到时间t_k的预测，Δt是从时刻t_k-1到t_k的时间间隔，

E (q) = [\begin{matrix} q_{4} & - q_{3} & q_{2} \\ q_{3} & q_{4} & - q_{1} \\ - q_{2} & q_{1} & q_{4} \\ - q_{1} & - q_{2} & - q_{3} \end{matrix}];

(3)读取星敏感器输出的光轴在惯性空间的指向

和利用步骤(2)更新的惯性姿态四元数将

和

从惯性坐标系转换到探测器本体系中，对应得到和

利用

和

和星敏感器的安装位置决定的星敏感器三轴指向

和

获取姿态预估的误差四元数

上标i表示惯性系，

Δ {\tilde{q}}_{k} = [\begin{matrix} δ \tilde{q} \\ 1 \end{matrix}] = [\begin{matrix} \frac{1}{2} (X_{star}^{b} \times {\hat{X}}_{star}^{b} + Y_{star}^{b} \times {\hat{Y}}_{star}^{b} + Z_{star}^{b} \times {\hat{Z}}_{star}^{b}) \\ 1 \end{matrix}];

(4)利用步骤(2)更新的惯性姿态四元数，计算得到陀螺漂移估计的残差

Δ {\tilde{b}}_{k}^{b} = {\tilde{ω}}_{k}^{b} - {\hat{b}}_{k | k - 1}^{b} - A ({\hat{q}}_{k | k - 1}) ω_{m}^{i},

其中b表示陀螺的常值漂移，

为月球的自转角速度矢量，A(q)表示由姿态四元数计算的姿态矩阵，

A (q) = [\begin{matrix} q_{1}^{2} - q_{2}^{2} - q_{3}^{2} + q_{4}^{2} & 2 (q_{1} q_{2} + q_{3} q_{4}) & 2 (q_{1} q_{3} - q_{2} q_{4}) \\ 2 (q_{1} q_{2} - q_{3} q_{4}) & - q_{1}^{2} + q_{2}^{2} - q_{3}^{2} + q_{4}^{2} & 2 (q_{2} q_{3} + q_{1} q_{4}) \\ 2 (q_{1} q_{3} + q_{2} q_{4}) & 2 (q_{2} q_{3} - q_{1} q_{4}) & - q_{1}^{2} - q_{2}^{2} + q_{3}^{2} + q_{4}^{2} \end{matrix}];

(5)采用卡尔曼滤波的方式获得星敏感器的姿态测量误差和陀螺常值漂移误差，其中状态方程和测量方程分别为：

[\begin{matrix} δ \overset{\cdot}{\hat{q}} \\ Δ {\dot{\hat{b}}}^{b} \end{matrix}] = [\begin{matrix} - [{\hat{ω}}^{b} \times] & - \frac{1}{2} I_{3 \times 3} \\ 0_{3 \times 3} & 0_{3 \times 3} \end{matrix}] [\begin{matrix} δ \hat{q} \\ Δ {\hat{b}}^{b} \end{matrix}] + [\begin{matrix} \frac{1}{2} I_{3 \times 3} & 0_{3 \times 3} \\ 0_{3 \times 3} & I_{3 \times 3} \end{matrix}] [\begin{matrix} η_{1}^{b} \\ η_{2}^{b} \end{matrix}]

[\begin{matrix} δ \tilde{q} \\ Δ {\tilde{b}}^{b} \end{matrix}] = [\begin{matrix} I_{3 \times 3} & 0_{3 \times 3} \\ 0_{3 \times 3} & I_{3 \times 3} \end{matrix}] [\begin{matrix} δ \hat{q} \\ Δ {\hat{b}}^{b} \end{matrix}] + [\begin{matrix} v_{star} \\ 0_{3 \times 1} \end{matrix}]

其中Δb^b表示陀螺漂移的误差，表示

的斜对称阵，

和

是陀螺噪声，V_star是星敏感器噪声，I表示单位阵，由此得到姿态和陀螺漂移估计的修正量：K为常系数卡尔曼滤波的稳态增益；

(6)利用步骤(5)得到的修正量对探测器的姿态四元数以及陀螺的常值漂移进行反馈修正，

当姿态估计稳定后进入下一步；

(7)利用步骤(1)和步骤(2)的方法，结合陀螺测量数据，对探测器进行姿态外推，得到探测器的惯性姿态四元数；

(8)利用加速度计的测量量

以及步骤(7)得到的惯性姿态四元数，外推探测器相对惯性系的速度，得到探测器相对惯性系的速度预测值以及加速度计的零偏预测值；

(9)根据探测器相对惯性系的速度预测值解算得到探测器相对月面的速度

{\hat{v}}_{am, k | k - 1}^{g} = [C_{m}^{g} C_{i}^{m} (t_{k})] \cdot [{\hat{v}}_{a, k | k - 1}^{i} - ω_{m}^{i} \times r_{a, k}^{i}],

其中

是月心固联坐标系到月理坐标系的转换矩阵，

表示由惯性系到月心固联坐标系的转换矩阵，

是探测器在惯性系下的位置矢量；

(10)以步骤(9)获得的探测器相对月面的速度作为测量量，采用卡尔曼滤波的方式获得惯导速度、姿态误差的估计值，其中状态方程和测量方程分别为：

[\begin{matrix} δ {\overset{\cdot}{v}}_{u} \\ δ {\overset{\cdot}{v}}_{e} \\ δ {\overset{\cdot}{v}}_{n} \\ {\overset{\cdot}{Φ}}_{u} \\ {\overset{\cdot}{Φ}}_{e} \\ {\overset{\cdot}{Φ}}_{n} \\ {\overset{\cdot}{&dtri;}}_{u} \end{matrix}] = [\begin{matrix} 0 & 2 ω_{m} \cos L & 0 & 0 & g_{n} & - g_{e} & 1 \\ - 2 ω_{m} \cos L & 0 & 2 ω_{m} \sin L & - g_{n} & 0 & g_{u} & 0 \\ 0 & - 2 ω_{m} \sin L & 0 & g_{e} & - g_{u} & 0 & 0 \\ 0 & 0 & 0 & 0 & ω_{m} \cos L & 0 & 0 \\ 0 & 0 & 0 & - ω_{m} \cos L & 0 & ω_{m} \sin L & 0 \\ 0 & 0 & 0 & 0 & - ω_{m} \cos L & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}] [\begin{matrix} δ v_{u} \\ δ v_{e} \\ δ v_{n} \\ Φ_{u} \\ Φ_{e} \\ Φ_{n} \\ {&dtri;}_{u} \end{matrix}]

[\begin{matrix} {\tilde{v}}_{u} \\ {\tilde{v}}_{e} \\ {\tilde{v}}_{n} \end{matrix}] = [\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \end{matrix}] [\begin{matrix} δ v_{u} \\ δ v_{e} \\ δ v_{n} \\ Φ_{u} \\ Φ_{e} \\ Φ_{n} \\ {&dtri;}_{u} \end{matrix}]

其中δV_u、δV_e、δV_n分别表示惯导解算的探测器对月速度在天向、东向和北向的误差；φ_u、φ_e和φ_n分别表示惯导在天向、东向和北向的姿态误差角；g_u、g_e和g_n则表示月理坐标系下天向、东向和北向的重力加速度；

是加速度零偏的天向投影；L是探测器在月表的纬度，且有

由此得到惯导误差的滤波修正量，K_IMU为常系数卡尔曼滤波的稳态增益，

[\begin{matrix} δ {\hat{v}}_{u} \\ δ {\hat{v}}_{e} \\ δ {\hat{v}}_{n} \\ {\hat{Φ}}_{u} \\ {\hat{Φ}}_{e} \\ {\hat{Φ}}_{n} \\ {\hat{&dtri;}}_{u} \end{matrix}] = K_{IMU} [\begin{matrix} - {\tilde{v}}_{u} \\ - {\tilde{v}}_{e} \\ - {\tilde{v}}_{n} \end{matrix}];

(11)利用步骤(10)获得的滤波估计值计算

由

解算出修正后的惯导姿态四元数

为探测器本体系相对月理坐标系的姿态矩阵，

{\hat{C}}_{g}^{g'} = [\begin{matrix} 1 & {\hat{Φ}}_{n} & - {\hat{Φ}}_{e} \\ - {\hat{Φ}}_{n} & 1 & {\hat{Φ}}_{u} \\ {\hat{Φ}}_{e} & - {\hat{Φ}}_{u} & 1 \end{matrix}];

(12)利用步骤(10)得到的探测器相对月面速度的修正量更新探测器的对月速度

{\hat{v}}_{am, k}^{g} = {\hat{v}}_{am, k | k - 1}^{g} + {[\begin{matrix} δ {\hat{v}}_{u} & δ {\hat{v}}_{e} & δ {\hat{v}}_{n} \end{matrix}]}^{T},

然后修正惯导外推的探测器惯性速度

并完成对加速度计零偏的预估值进行修正，

{\hat{&dtri;}}_{k - 1}^{b} = {\hat{&dtri;}}_{k | k - 1}^{b} + {\hat{C}}_{g}^{b} (t_{k}) {\hat{C}}_{g}^{g^{'}} [\begin{matrix} {\hat{&dtri;}}_{u} \\ 0 \\ 0 \end{matrix}] .

本发明与现有技术相比的优点在于：

(1)本发明方法在对准过程中使用通用和自动的星敏感器取代了人工且复杂的瞄准设备，使得月面惯导对准可以在无人干预的情况下实施，简化了系统构成，具有成本低、适用面广、自主性强的特点；同时，也解决了月球自转角速度低造成惯导自对准无法达到足够的方位对准精度的问题。

(2)本发明方法改变了常规将各类数据同时引入修正的思路，将方法分成星敏感器+陀螺姿态滤波和惯导误差修正两个顺序进行的步骤。在第一步中用星敏感器和陀螺确定探测器姿态的初值，并标定陀螺漂移；在第二步中校正由星敏和IMU之间安装误差引起的惯导平台误差，并估计天向加速度计零偏。这一操作降低了滤波方程的维数，简化了计算量，便于星上实施，同时也解决了起飞前最后阶段星敏感器无法使用(为了防止发动机吹起月尘污染镜头，星敏感器在起飞十几分钟前需要关闭防尘盖，因此步骤一无法进行到探测器点火起飞前)的问题。

附图说明

图1为本发明方法的原理框图；

图2为本发明方法中月面起飞对准姿态误差角变化曲线。

具体实施方式

本发明方法在具体实施上分为两个主要的步骤，具体如图1所示。

1、惯性姿态估计

使用星敏感器+陀螺进行惯性姿态估计，即首先使用陀螺测量预估停留在月面的探测器惯性姿态，然后利用星敏感器测量对陀螺预估的姿态进行修正，同时估计出陀螺角速度测量误差。该估计方法与在轨飞行器常规的星敏感器+陀螺定姿方法的不同之处为：当探测器停留在月面时，它的角速度等于月球自转角速度，这个量已知，因此月球自转角速度也可以用作陀螺漂移估计的参考。星敏感器+陀螺获得的惯性姿态在滤波稳定后赋给惯导系统建立初值。该步骤的具体操作如下：

a)姿态外推

假设当前时刻t_k陀螺获得探测器的角速度为而前一次的陀螺漂移估计为

由于陀螺常值漂移不随时间变化，则t_k时刻探测器的角速度的估计值

为

{\hat{ω}}_{k}^{b} = {\tilde{ω}}_{k}^{b} - {\hat{b}}_{k | k - 1}^{b} - - - (1)

其中，

{\hat{b}}_{k | k - 1}^{b} = {\hat{b}}_{k - 1}^{b} - - - (2)

以上各式中，上标b表示该矢量描述在探测器本体坐标系中，下标k表示时间t_k，k|k-1表示由时间t_k-1到时间t_k的预测，波浪号表示测量值，上三角表示估计值，下同。

由此，可以根据姿态运动学方程，对探测器的惯性姿态进行预估。姿态运动学方程为

\overset{\cdot}{q} = \frac{1}{2} E (q) ω - - - (3)

其中

E (q) = [\begin{matrix} q_{4} & - q_{3} & q_{2} \\ q_{3} & q_{4} & - q_{1} \\ - q_{2} & q_{1} & q_{4} \\ - q_{1} & - q_{2} & - q_{3} \end{matrix}] - - - (4)

q是表征探测器本体系相对惯性系的姿态四元数，且有q＝[q₁ q₂ q₃ q₄]^T。q上方的一点表示q对时间的导数，ω为探测器的角速度。

对(3)式进行离散化，就可以使用如下方程完成从时刻t_k-1到时刻t_k的姿态四元数的预估

{\hat{q}}_{k | k - 1} \approx {\hat{q}}_{k - 1} + \frac{1}{2} E ({\hat{q}}_{k - 1}) {\hat{ω}}_{k}^{b} \cdot Δt - - - (5)

Δt是从时刻t_k-1到t_k的时间间隔，即Δt＝t_k-t_k-1。

b)常系数滤波1

星敏感器能够输出其敏感器光轴和横轴在惯性空间的指向，即

和(上标i表示该矢量描述在惯性系下，下标star表示星敏感器)，利用预估的姿态四元数可以将它们转换到探测器的本体系中，

\{\begin{matrix} {\hat{X}}_{star}^{b} = A ({\hat{q}}_{k | k - 1}) {\tilde{X}}_{star}^{i} \\ {\hat{Y}}_{star}^{b} = A ({\hat{q}}_{k | k - 1}) {\tilde{Y}}_{star}^{i} \\ {\hat{Z}}_{star}^{b} = A ({\hat{q}}_{k | k - 1}) {\hat{Z}}_{star}^{i} \end{matrix} - - - (6)

其中A(q)表示将四元数q转换为姿态矩阵

A (q) = [\begin{matrix} q_{1}^{2} - q_{2}^{2} - q_{3}^{2} + q_{4}^{2} & 2 (q_{1} q_{2} + q_{3} q_{4}) & 2 (q_{1} q_{3} - q_{2} q_{4}) \\ 2 (q_{1} q_{2} - q_{3} q_{4}) & - q_{1}^{2} + q_{2}^{2} - q_{3}^{2} + q_{4}^{2} & 2 (q_{2} q_{3} + q_{1} q_{4}) \\ 2 (q_{1} q_{3} + q_{2} q_{4}) & 2 (q_{2} q_{3} - q_{1} q_{4}) & - q_{1}^{2} - q_{2}^{2} + q_{3}^{2} + q_{4}^{2} \end{matrix}] - - - (7)

由于星敏感器的安装已知，即确切的星敏感器三轴指向在探测器本体系是已知的，它们分别为

和

那么根据这三个轴的理论值与预估值的差可以获得姿态预估误差四元数的测量值

Δ {\tilde{q}}_{k} \approx [\begin{matrix} δ \tilde{q} \\ 1 \end{matrix}] = [\begin{matrix} \frac{1}{2} (X_{star}^{b} \times {\hat{X}}_{star}^{b} + Y_{star}^{b} \times {\hat{Y}}_{star}^{b} + Z_{star}^{b} \times {\hat{Z}}_{star}^{b}) \\ 1 \end{matrix}] - - - (8)

由于月球自转角速度是已知的，记为

那么也可以根据预估的姿态四元数计算出陀螺漂移估计的残差

Δ {\tilde{b}}_{k}^{b} = {\tilde{ω}}_{k}^{b} - {\hat{b}}_{k | k - 1}^{b} - A ({\hat{q}}_{k | k - 1}) ω_{m}^{i} - - - (9)

由于姿态预估误差四元数为

对其两边求导并将式(3)代入其中，则可以推导出姿态四元数外推过程中标量部分的误差传播方程为

δ \overset{\cdot}{q} = - {\hat{ω}}^{b} \times δq - \frac{1}{2} Δ b^{b} - \frac{1}{2} η_{1}^{b} - - - (10)

式中

是陀螺测量的随机白噪声。而陀螺的常值漂移估计误差可以建模为

Δ {\overset{\cdot}{b}}^{b} = η_{2}^{b} - - - (11)

也是随机白噪声。联立式(10)、(11)可以构造滤波器的状态方程为

[\begin{matrix} δ \overset{\cdot}{\hat{q}} \\ Δ {\dot{\hat{b}}}^{b} \end{matrix}] = [\begin{matrix} - [{\hat{ω}}^{b} \times] & - \frac{1}{2} I_{3 \times 3} \\ 0_{3 \times 3} & 0_{3 \times 3} \end{matrix}] [\begin{matrix} δ \hat{q} \\ Δ {\hat{b}}^{b} \end{matrix}] + [\begin{matrix} \frac{1}{2} I_{3 \times 3} & 0_{3 \times 3} \\ 0_{3 \times 3} & I_{3 \times 3} \end{matrix}] [\begin{matrix} η_{1}^{b} \\ η_{2}^{b} \end{matrix}] - - - (12)

其中

表示

的斜对称阵，即

[{\hat{ω}}^{b} \times] = [\begin{matrix} 0 & - {\hat{ω}}_{z}^{b} & {\hat{ω}}_{y}^{b} \\ {\hat{ω}}_{z}^{b} & 0 & - {\hat{ω}}_{x}^{b} \\ - {\hat{ω}}_{y}^{b} & {\hat{ω}}_{x}^{b} & 0 \end{matrix}]

(下标xyz分别表示

在探测器本体系三个坐标轴的分量，探测器本体系固联在探测器上，原点在探测器质心，三个坐标轴指向探测器上固定的方向。)

滤波器的测量方程为

[\begin{matrix} δ \tilde{q} \\ Δ {\tilde{b}}^{b} \end{matrix}] = [\begin{matrix} I_{3 \times 3} & 0_{3 \times 3} \\ 0_{3 \times 3} & I_{3 \times 3} \end{matrix}] [\begin{matrix} δ \hat{q} \\ Δ {\hat{b}}^{b} \end{matrix}] + [\begin{matrix} v_{star} \\ 0_{3 \times 1} \end{matrix}] - - - (13)

其中，V_star表示星敏感器的测量噪声，I_3×3是3×3的单位矩阵，0_3×3是3×3的零矩阵。式(13)中的测量量由式(8)和(9)得到。由此，可以使用卡尔曼滤波估计出姿态误差和陀螺常值漂移误差。

令系统状态

x = [\begin{matrix} δq \\ Δ b^{b} \end{matrix}] - - - (14)

根据式(12)状态转移矩阵为

Φ_{k | k - 1} \approx I + Δt [\begin{matrix} - [{\hat{ω}}_{k - 1}^{b} \times] & - \frac{1}{2} I_{3 \times 3} \\ 0_{3 \times 3} & 0_{3 \times 3} \end{matrix}] - - - (15)

根据式(13)观测矩阵为

H = [\begin{matrix} I_{3 \times 3} & 0_{3 \times 3} \\ 0_{3 \times 3} & I_{3 \times 3} \end{matrix}] - - - (16)

系统噪声

w = [\begin{matrix} η_{1}^{b} \\ η_{2}^{b} \end{matrix}] - - - (17)

其方差为Q。并且令

G = [\begin{matrix} \frac{1}{2} I_{3 \times 3} & 0_{3 \times 3} \\ 0_{3 \times 3} & I_{3 \times 3} \end{matrix}] - - - (18)

令测量量为

z = [\begin{matrix} δ \tilde{q} \\ Δ {\tilde{b}}^{b} \end{matrix}]

测量噪声为

v = [\begin{matrix} v_{star} \\ 0_{3 \times 1} \end{matrix}] - - - (19)

它的方差为R，则可以使用卡尔曼滤波进行状态估计，方法如下

x_k|k-1＝Φ_k|k-1x_k-1 (20)

P_{k | k - 1} = Φ_{k | k - 1} P_{k - 1} Φ_{k | k - 1}^{T} + G_{k} Q_{k} G_{k}^{T} - - - (21)

K_{k} = {(P_{k | k - 1} H_{k}^{T})}^{- 1} (H_{k} P_{k | k - 1} H_{k}^{T} + R_{k}) - - - (22)

x_k＝x_k|k-1+K_k(z_k-H_kx_k|k-1) (23)

P_{k} = (I - K_{k} H_{k}) P_{k | k - 1} {(I - K_{k} H_{k})}^{T} + K_{k} R_{k} K_{k}^{T} - - - (24)

从式(20)～(24)看到矩阵K_k的稳态值(k→∞)只与Φ、G、Q、H和R矩阵有关，除了Φ以外，都是常数矩阵；而且由于探测器停留在月面时自转角速度不变，因此用月球自转角速度矢量在标称姿态下的投影

取代

后，Φ也是常数矩阵。则该滤波可以转换为常系数滤波，即

首先，用式(8)和(9)式获得观测残差

和

然后，由观测残差获得姿态和陀螺漂移估计的修正量

[\begin{matrix} δ {\hat{q}}_{k} \\ Δ {\hat{b}}_{k}^{b} \end{matrix}] = K [\begin{matrix} δ {\tilde{q}}_{k} \\ Δ {\tilde{b}}_{k}^{b} \end{matrix}] - - - (25)

其中滤波修正系数矩阵K取值为由式(21)、(22)和(24)递推得到的卡尔曼滤波稳态增益，即k足够大时的K_k。

最后，进行闭环修正

{\hat{q}}_{k} = {\hat{q}}_{k | k - 1} + E (\hat{q}) δ {\hat{q}}_{k} - - - (26)

{\hat{b}}_{k}^{b} = {\hat{b}}_{k | k - 1}^{b} + Δ {\hat{b}}_{k}^{b} - - - (27)

2、惯性平台误差的估计和校正

在惯性姿态估计完成后，启动惯性导航，利用加速度计和陀螺测量进行探测器速度、姿态外推。由于存在IMU测量误差、星敏感器测量误差以及星敏感器和IMU之间的安装误差，惯导递推获得的探测器速度和姿态也存在误差并且误差逐渐发散。但由于探测器相对月表并没有运动，因此可以将惯导解算出的对月速度作为观测量，估计惯导平台的误差角以及陀螺和加速度计的测量误差，即自对准。但在这一过程中，滚动误差角与东向陀螺漂移、俯仰误差角与北向加速度计零偏、偏航误差角与东向加速度计零偏这三对组合中的两个状态量不同时可观，因此滤波状态中需要去掉这六个状态中的三个。另外，由于陀螺漂移已经进行了估计，不需在在此进行重复估计。因此，最终滤波状态中去掉了三个陀螺漂移和北向及东向加速度计零偏。

由滤波器估计出惯导平台的误差角后，将该误差角反馈给惯导系统，即可校正惯性导航外推的探测器速度和姿态。

a)姿态外推

步骤二的姿态外推与步骤一的姿态外推基本相同，区别是陀螺漂移估计不再随时间更新。即

{\hat{ω}}_{k}^{b} = {\tilde{ω}}_{k}^{b} - {\hat{b}}^{b} - - - (28)

{\hat{q}}_{k | k - 1} \approx {\hat{q}}_{k - 1} + \frac{1}{2} E ({\hat{q}}_{k - 1}) {\hat{ω}}_{k}^{b} \cdot Δt - - - (29)

其中，是步骤一最终获得的陀螺常值漂移估计。

b)速度外推

取参考坐标系为月心惯性坐标系，探测器的速度矢量表示为

假设加速度计的测量为

加速度计零偏的估计值为

引力加速度为g，则从时间t_k-1到t_k的速度外推可如下计算：

{\hat{v}}_{a, k | k - 1}^{i} = {\hat{v}}_{a, k - 1}^{i} + Δt [A ({\hat{q}}_{k - 1}) \cdot ({\hat{f}}_{k - 1}^{b} - {\hat{&dtri;}}_{k - 1}^{b}) - g_{k - 1}^{i}] - - - (30)

同时对加速度计零偏进行预测

{\hat{&dtri;}}_{k | k - 1}^{b} = {\hat{&dtri;}}_{k - 1}^{b} - - - (31)

c)预估探测器相对月面的速度

首先将每一时刻陀螺预估的姿态四元数转换为惯性姿态矩阵

{\hat{C}}_{i}^{b} (t_{k}) = A ({\hat{q}}_{k}) - - - (32)

由于月球的星历已知，因此可以得到从惯性坐标系到月心固联坐标系的变换阵

月心固联坐标系的原点在月心，X轴在赤道面内指向零度经线，Z轴平行于月球自转轴指向北极，Y轴在赤道面内与X轴、Y轴构成右手坐标系。由于探测器在月心固联系下的位置已知，设经度、纬度和高度分别为L、λ和h，由此可以计算出月心固联系到当地月理坐标系(原点在当地参考椭球面上，X轴垂直于当地参考椭球面指向天，Y轴指向东，Z轴指向北)的转换阵。

C_{m}^{g} = [\begin{matrix} \cos L \cos λ & \cos L \sin λ & \sin L \\ - \sin & \cos λ & 0 \\ - \sin L \cos λ & - \sin L \sin λ & \cos L \end{matrix}] - - - (33)

这样，就可以算出探测器本体相对当地月理坐标系的姿态矩阵

{\hat{C}}_{g}^{b} (t_{k}) = {\hat{C}}_{i}^{b} (t_{k}) {[C_{m}^{g} C_{i}^{m} (t_{k})]}^{T} - - - (34)

然后，根据惯性速度解算探测器相对月面的速度

{\hat{v}}_{am, k | k - 1}^{g} = [C_{m}^{g} C_{i}^{m} (t_{k})] \cdot [{\hat{v}}_{a, k | k - 1}^{i} - ω_{m}^{i} \times r_{a, k}^{i}] - - - (35)

其中，

是月球自转角速度矢量，是探测器在惯性系下的位置，它可以如下计算

r_{a, k}^{i} = {[C_{m}^{g} C_{i}^{m} (t_{k})]}^{T} [\begin{matrix} R_{m} + h \\ 0 \\ 0 \end{matrix}] - - - (36)

其中，R_m为月球参考半径。

d)常系数滤波2

由于探测器停留在月面不动，真实的相对月面的速度为零，因此惯导系统外推得到的探测器相对月面的速度就是速度误差。这样可建立如下滤波系统方程和测量方程

[\begin{matrix} δ {\overset{\cdot}{v}}_{u} \\ δ {\overset{\cdot}{v}}_{e} \\ δ {\overset{\cdot}{v}}_{n} \\ {\overset{\cdot}{Φ}}_{u} \\ {\overset{\cdot}{Φ}}_{e} \\ {\overset{\cdot}{Φ}}_{n} \\ {\overset{\cdot}{&dtri;}}_{u} \end{matrix}] = [\begin{matrix} 0 & 2 ω_{m} \cos L & 0 & 0 & g_{n} & - g_{e} & 1 \\ - 2 ω_{m} \cos L & 0 & 2 ω_{m} \sin L & - g_{n} & 0 & g_{u} & 0 \\ 0 & - 2 ω_{m} \sin L & 0 & g_{e} & - g_{u} & 0 & 0 \\ 0 & 0 & 0 & 0 & ω_{m} \cos L & 0 & 0 \\ 0 & 0 & 0 & - ω_{m} \cos L & 0 & ω_{m} \sin L & 0 \\ 0 & 0 & 0 & 0 & - ω_{m} \cos L & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}] [\begin{matrix} δ v_{u} \\ δ v_{e} \\ δ v_{n} \\ Φ_{u} \\ Φ_{e} \\ Φ_{n} \\ {&dtri;}_{u} \end{matrix}] - - - (37)

[\begin{matrix} {\tilde{v}}_{u} \\ {\tilde{v}}_{e} \\ {\tilde{v}}_{n} \end{matrix}] = [\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \end{matrix}] [\begin{matrix} δ v_{u} \\ δ v_{e} \\ δ v_{n} \\ Φ_{u} \\ Φ_{e} \\ Φ_{n} \\ {&dtri;}_{u} \end{matrix}]

其中，δV_u、δV_e、δV_n分别表示惯导解算的探测器对月速度在天向、东向和北向的误差；φ_u、φ_e和φ_n分别表示惯导在天向、东向和北向的姿态误差角；g_u、g_e和g_n则表示月理坐标系下天向、东向和北向的重力加速度；

是加速度零偏的天向投影。由于探测器相对月表没有运动，因此惯导系统给出的对月速度的误差就是惯导给出的对月速度本身，因此对于t_k时刻的测量有

{[\begin{matrix} {\tilde{v}}_{u, k} & {\tilde{v}}_{e, k} & {\tilde{v}}_{n, k} \end{matrix}]}^{T} = {\hat{v}}_{am, k | k - 1}^{g} - - - (38)

式(37)表示的是一个定常系统，因此可以使用常系数滤波进行状态估计。

令

δX = [\begin{matrix} {δv}^{g} & Φ \end{matrix}] = {[\begin{matrix} {δv}_{u} & {δv}_{e} & {δv}_{n} & Φ_{u} & Φ_{e} & Φ_{n} & {&dtri;}_{u} \end{matrix}]}^{T} - - - (39)

则滤波方程可以转换为

δ \hat{X} = K_{IMU} [\begin{matrix} - {\tilde{v}}_{u} \\ - {\tilde{v}}_{e} \\ - {\tilde{v}}_{n} \end{matrix}] - - - (40)

其中，K_IMU是式(37)构成的卡尔曼滤波的稳态增益阵，它可以用与式(21)、(22)和(24)同样的方法获得。

e)惯性姿态修正

获得月理系下的姿态误差后，需要将它转换为惯性系下的姿态误差，并反馈修正惯导姿态四元数。修正后的惯性姿态矩阵为

{\hat{C}}_{i}^{b} (t_{k}) = {\hat{C}}_{g}^{b} (t_{k}) {\hat{C}}_{g}^{g'} {\hat{C}}_{m}^{g} C_{i}^{m} (t_{k}) - - - (41)

其中

{\hat{C}}_{g}^{g'} = [\begin{matrix} 1 & {\hat{Φ}}_{n} & - {\hat{Φ}}_{e} \\ - {\hat{Φ}}_{n} & 1 & {\hat{Φ}}_{u} \\ {\hat{Φ}}_{e} & - {\hat{Φ}}_{u} & 1 \end{matrix}] - - - (42)

由

可以解算出修正后的姿态四元数

f)惯性速度修正

首先获得相对月面速度的修正量

{\hat{v}}_{am, k}^{g} = {\hat{v}}_{am, k | k - 1}^{g} + {[\begin{matrix} δ {\hat{v}}_{u} & δ {\hat{v}}_{e} & δ {\hat{v}}_{w} \end{matrix}]}^{T} - - - (43)

然后，修正探测器惯性速度

{\hat{v}}_{a, k}^{i} = {[C_{m}^{g} C_{i}^{m} (t_{k})]}^{T} {\hat{v}}_{am, k}^{g} + ω_{m}^{i} \times r_{a, k}^{i} - - - (44)

之后，还要对加速度计零偏的预估值进行修正

{\hat{&dtri;}}_{k - 1}^{b} = {\hat{&dtri;}}_{k | k - 1}^{b} + {\hat{C}}_{g}^{b} (t_{k}) {\hat{C}}_{g}^{g^{'}} [\begin{matrix} {\hat{&dtri;}}_{u} \\ 0 \\ 0 \end{matrix}] - - - (45)

实施例

此处通过一个仿真例子对本发明提出的月面起飞对准方法进行验证。设探测器在月面的位置为北纬45°、西经60°、高度0m，且探测器纵轴朝上，但相对月面有15°的倾斜。地面给定的探测器位置在经度、纬度和高度上均有500m(3σ)的误差。

探测器上装备有星敏感器和IMU。其中，星敏感器光轴的标称指向为[cos(38°)，cos(54．7°)，cos(102．2°)]^T，光轴误差为3"、横轴误差为24"(3σ)；IMU包含三个陀螺和三个加速度计，它们的输入轴指向相同。分别为

P₁[cos(54°44′8")，sin(54°44′8")cos(90°)，sin(54°44′8")sin(90°)]^T

P₂=[cos(54°44′8")，sin(54°44′8")cos(210°)，sin(54°44′8")sin(210°)]^T

P₃=[cos(54°44′8")，sin(54°44′8″)cos(330°)，sin(54°44′8")sin(330°)]^T陀螺的常值漂移为3°/h(3σ)，加速度计的零偏为1×10^-4g(3σ)。

考虑到探测器经过火箭发射、月球软着陆等飞行阶段后，内部结构存在变形，假设月面起飞前星敏相对IMU在光轴和横轴方向分别有20"和2′的结构变形。

仿真中设定对准总共用时1200秒，其中前300秒进行步骤一，后900秒进行步骤二。对准的姿态误差变化曲线如图2所示。可以看到，在步骤一中，通过星敏感器和陀螺滤波估计，惯导姿态能够很快的稳定。但由于星敏感器和IMU之间存在的形变误差，使得第一步的对准结果存在常值偏差。到步骤二后，由于引入了加速度计测量作为重力方向的参考，因此能够校正星敏感器和IMU安装偏差在水平方向带来的对准误差，即俯仰θ和偏航ψ误差角明显减小。仿真表明，该月面对准方法简单、实用且有效。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。