CN105737858A - 一种机载惯导系统姿态参数校准方法与装置 - Google Patents

Abstract

本发明公开了一种机载惯导系统姿态参数校准方法与装置，以高动态星敏感器为数据获取手段，以恒星为测量对象，通过探测天球上的多颗恒星，并进行星点定位、星图识别和姿态解算，提供星敏感器相对于惯性参考系的姿态；然后，根据星敏感器坐标系与惯导系统坐标系间的变换关系，得到惯导坐标系下的飞机姿态基准数据；最后，以此基准数据为参考，与惯导系统的姿态数据进行比对，从而实现机载环境下对惯导系统姿态精度的动态在线校准。

Description

一种机载惯导系统姿态参数校准方法与装置

技术领域

本发明属于导航姿态参数校准的技术领域，具体涉及一种机载惯导系统姿态参数校准方法与装置。

背景技术

惯性导航系统是一种不依赖于外部信息、也不向外部辐射能量的自主式导航系统，具有很好的隐蔽性，不受外界电磁干扰的影响。在航空领域中，它是各种飞机不可缺少的导航系统，其姿态精度性能对飞机导航至关重要，因此必须对惯导系统的姿态精度进行校准。与其它机载设备一样，惯性导航系统必须在真实飞行条件下进行试飞，从而为新型惯性导航系统的设计定型/鉴定提供技术依据，但是，长期以来，惯性导航系统的试飞鉴定、特别是实际使用环境中姿态参数的校准问题尚未很好解决。

目前国内一般采用地面转台对惯导系统参数进行校准，中国专利“一种POS方位精度和姿态精度的地面测试方法”(CN201010613323.X)和中国专利“一种高精度多惯导系统姿态精度评定方法”(CN201210156278.9)均是在地面转台上进行，通过设定轨迹并旋转转台，以转台提供的方位角和水平姿态角信息为基准，对惯导系统输出姿态参数进行精度校准。地面转台试验可以对惯导系统进行高精度校准，但是，实际中飞机受到天气、气流等影响，导致惯导系统工作在复杂的动态环境中，这种使用环境与实验室校准平台环境存在较大差异，从而对实验室校准结果的可信度产生影响。

星敏感器以天球惯性坐标系中的恒星为参考基准，输出载体在惯性空间中的绝对姿态信息，具有无漂移、精度高等特点，但是，当用于飞机机载设备时，飞机机动性较强，传统星敏感器由于曝光时间较长，出现成像光斑拖尾的现象，造成星点定位精度下降，严重时甚至不能输出姿态。相比于传统星敏感器，高动态星敏感器具有极高的探测灵敏度，曝光时间相对于传统星敏感器大为缩短，在实际飞行的复杂动态环境下仍然能够实现高精度姿态输出，因此本发明以高动态星敏感器输出姿态数据为基准，完成机载惯导系统姿态参数的动态在线校准。

发明内容

本发明提出一种基于高动态星敏感器的惯导系统姿态参数的校准方法与装置，以高动态星敏感器为数据获取手段，以恒星为测量对象，通过探测天球上的多颗恒星，并进行星点定位、星图识别和姿态解算，提供星敏感器相对于惯性参考系的姿态；然后，根据星敏感器坐标系与惯导系统坐标系间的变换关系，得到惯导坐标系下的飞机姿态基准数据；最后，以此基准数据为参考，与惯导系统的姿态数据进行比对，从而实现机载环境下对惯导系统姿态精度的动态在线校准。

本发明采用的技术方案为：一种机载惯导系统姿态参数校准装置，该校准装置包括高动态星敏感器、GPS同步授时系统、数据采集器、稳压电源和固定支架，其中，固定支架用于固定安装校准装置高动态星敏感器、GPS同步授时系统、数据采集器和稳压电源以及待校准的惯导系统，稳压电源则为高动态星敏感器、GPS同步授时系统和数据采集器以及待校准的惯导系统供电，当需要校准时，将待校准的惯导系统安装在固定支架上，高动态星敏感器将姿态数据及其对应时刻的脉冲信号输出到数据采集器，GPS同步授时系统用于提供标准GPS时间，其中，待校准惯导系统数据直接进行授时，而高动态星敏感器则需要首先对脉冲信号进行授时，然后利用脉冲信号与高动态星敏感器姿态数据一一对应的性质，间接完成对高动态星敏感器输出姿态数据的授时。

本发明还提供一种机载惯导系统姿态参数校准方法，利用上述的机载惯导系统姿态参数校准装置，该校准方法步骤如下：

步骤一、坐标系基准统一步骤；

将星敏感器与惯导系统的相对参考系统一，即建立惯性坐标系与地理坐标系的转换关系；然后，利用惯导系统姿态数据与星敏感器姿态数据将二者间的常值误差矩阵解算出来；

步骤二、数据采集与校准步骤；

步骤(1)、以恒星为测量对象，通过高动态星敏感器探测天球上的多颗恒星，并进行星点定位、星图识别和姿态解算，提供星敏感器相对于惯性参考系的姿态；以步骤一得到的星敏感器坐标系与惯导系统坐标系的基准统一变换关系，得到惯导坐标系下的飞机姿态基准数据；

步骤(2)、以步骤(1)得到的基准数据为参考，与惯导系统的姿态数据进行比对，从而实现机载环境下对惯导系统姿态精度的动态在线校准。

其中，步骤一坐标系基准统一的关键在于常值误差矩阵B_s ^g的解算，通过选取N个采样点，每个采样点采集M对有效数据，则共有N×M个星敏感器与惯导系统的有效数据对，每对数据满足式(3)，式中，表示第i个采样点的第j组惯导系统相对于地理坐标系的姿态数据与星敏感器相对于地理坐标系的姿态数据组成的数据对，B_s ^g为待求的常值误差矩阵，

$Q_{t\left(ij\right)}^g=B_s^g{Q}_{t\left(ij\right)}^s,1\≤i\≤N,1\≤j\≤M---\left(3\right)$

然后，采用QUEST算法求解B_s ^g的最佳估计值，其步骤如下：

①首先，对式(3)进行改写，其等价表达式如式(4)所示，式中，r₁、r₂、r₃和w₁、w₂、w₃分别表示的列向量，

${\left[\begin{array}{ccc}{r}_1& r_2& r_3\end{array}\right]}_{\left(ij\right)}=B_s^g{\left[\begin{array}{ccc}{w}_1& w_2& w_3\end{array}\right]}_{\left(ij\right)},1\≤i\≤N,1\≤j\≤M---\left(4\right)$

进一步的，将式(4)写成列向量分量形式如式(5)所示，式中(r_ijk,w_ijk)表示第i个采样点的第j组对应姿态数据的第k个列向量分量对，

$r_{ijk}=B_s^g{w}_{ijk},1\≤i\≤N,1\≤j\≤M,1\≤k\≤3---\left(5\right)$

②经过上述改写后，式(3)等价为式(5)，而式(5)为典型的Wahba问题，即求解B_s ^g的最佳估计值，等效为求解行列式为1的最优正交矩阵B_s ^g，使得损失函数：

$L\left(B_s^g\right)\≡\frac{1}{2}\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}{\left|r_{ijk}-B_s^g{w}_{ijk}\right|}^2---\left(6\right)$

最小，式中，a_ijk为非负系数，这里取a_ijk恒为1；

③最后，利用QUEST算法，求解B_s ^g的最优估计值：

首先，由获得的有效数据，分别计算有关中间变量，其公式如式(7)～(11)所示。

$B=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}{r}_{ijk}{\left(w_{ijk}\right)}^T---\left(7\right)$

$\mathrm{\σ}=tr\left(B\right)=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}{r}_{ijk}\·w_{ijk}---\left(8\right)$

S＝B+B^T(9)

$Z=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}(r_{ijk}\×w_{ijk})---\left(10\right)$

利用上述中间变量，可以将矩阵B_s ^g的最优估计问题，转化为与之等价的四元数的最优估计问题，其公式如(12)所示，式中λ_max为矩阵K的最大特征值，与B_s ^g等价，为λ_max对应的特征向量，

$K{\stackrel{\‾}{q}}_{opt}={\mathrm{\λ}}_{max}{\stackrel{\‾}{q}}_{opt}---\left(12\right)$

对式(12)，在给定精度下可以利用牛顿-拉夫逊迭代法求解则由及式(7)～(10)可得，的最优估计计算过程如式(13)～(16)所示，其中，式(13)中的tr(A)、adj(A)和det(A)分别表示矩阵A的迹、伴随矩阵和行列式，

κ＝tr(adj(S)),△＝det(S)(13)

$\left\{\begin{array}{c}\mathrm{\α}={\left({\widehat{\mathrm{\λ}}}_{max}\right)}^2-{\mathrm{\σ}}^2+\mathrm{\κ}\\ \mathrm{\β}={\widehat{\mathrm{\λ}}}_{max}-\mathrm{\σ}\\ \mathrm{\γ}=({\widehat{\mathrm{\λ}}}_{max}+\mathrm{\σ})\·\mathrm{\α}-\mathrm{\Δ}\end{array}\right.---\left(14\right)$

X＝(αI+βS+S²)Z(15)

${\stackrel{\‾}{q}}_{opt}={\stackrel{\→}{q}}_r+q_0={(q_1,q_2,q_3,q_0)}^T=\frac{1}{\sqrt{{\mathrm{\γ}}^2+{\left|X\right|}^2}}\left\{\begin{array}{c}X\\ \mathrm{\γ}\end{array}\right\}---\left(16\right)$

式(16)所示就是的最优估计，最后，通过四元数与旋转矩阵的转换关系，可以得到B_s ^g的最优估计结果如式(17)所示，

${\hat{B}}_s^g=\left[\begin{array}{ccc}{q_0}^2+{q_1}^2-{q_2}^2-{q_3}^2& 2(q_1{q}_2+q_0{q}_3)& 2(q_1{q}_3-q_0{q}_2)\\ 2(q_1{q}_2-q_0{q}_3)& {q_0}^2-{q_1}^2+{q_2}^2-{q_3}^2& 2(q_3{q}_2+q_0{q}_1)\\ 2(q_1{q}_3+q_0{q}_2)& 2(q_3{q}_2-q_0{q}_1)& {q_0}^2-{q_1}^2-{q_2}^2+{q_3}^2\end{array}\right]---\left(17\right)$

其中，步骤二数据采集与校准，由于测量原理不同，现有星敏感器的姿态数据更新率一般低于惯导系统的更新率，且由于时间同步误差，星敏感器与惯导系统的姿态测量时刻不能完全一致，故而无法直接应用星敏感器数据对惯导系统姿态参数进行校准，而是需要利用星敏感器所获得的离散姿态真值数据点插值获得姿态真值曲线，然后利用真值曲线对惯导系统姿态参数进行校准与评定。

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

(1)本发明充分利用高动态星敏感器在复杂动态环境下仍然能够实现高精度姿态输出的特点，以高动态星敏感器输出姿态数据为基准，与惯导系统的姿态数据进行比对，从而实现机载环境下对惯导系统姿态精度的动态在线校准，提高了机载惯导系统姿态参数校准的可靠性与可信度；

(2)本发明通过选取N×M个星敏感器与惯导系统的有效数据对，利用QUEST算法对星敏感器与惯导系统之间的常值安装误差矩阵B_s ^g给出最优估计值，从而提高了惯导系统姿态参数校准的精度；

(3)本发明采用插值方法逼近真值曲线，完成了对惯导系统姿态参数的误差大小及数据偏差情况的校准与评定，克服了星敏感器的基准姿态数据与待校准惯导系统的姿态数据在时刻上不完全一致的问题。

附图说明

图1为惯导系统姿态参数校准装置示意图；

图2为机载惯导系统姿态参数校准方法流程示意图；

图3为各坐标系相互位置关系示意图；

图4为星敏感器与惯导系统坐标系误差关系示意图；

图5为姿态参数校准实验过程示意图。

具体实施方式

下面结合附图以及具体实施方式进一步说明本发明。

图1所示是本发明提出的机载惯导系统姿态参数校准装置示意图。由图1可知，整个校准装置由高动态星敏感器、GPS同步授时系统、数据采集器、稳压电源和固定支架等组成，其中，固定支架用于固定安装校准装置各部件，稳压电源则为各部件供电。当需要校准时，将待校准惯导系统安装在固定支架上，高动态星敏感器将姿态数据及其对应时刻的脉冲信号输出到数据采集器，GPS同步授时系统用于提供标准GPS时间，其中，待校准惯导系统数据直接进行授时，而高动态星敏感器则需要首先对脉冲信号进行授时，然后利用脉冲信号与高动态星敏感器姿态数据一一对应的性质，间接完成对高动态星敏感器输出姿态数据的授时。

以图1所示校准装置为基础，机载惯导系统的姿态参数校准方法主要分为两个步骤，如图2所示，即坐标系基准统一阶段和数据采集与校准阶段两步：

一、坐标系基准统一

(1)各坐标系说明

各坐标系的相互关系如图3所示，为了便于描述，对各坐标系说明如下：

1)地心赤道惯性坐标系O_eX_iY_iZ_i：该坐标系是惯性坐标系，其原点选取在地心，O_eZ_i轴指向地轴方向，O_eX_i轴与O_eY_i轴位于赤道平面内，O_eX_i轴指向春分点，且三轴构成右手坐标系，其中将2000年1月1日12时对应春分点所建立坐标系称为J2000.0惯性坐标系；

2)地球坐标系O_eX_eY_eZ_e：该坐标系与地球固连，因而始终相对于地球静止，其原点为地心，O_eZ_e轴指向地轴方向，O_eX_e轴与O_eY_e轴位于赤道平面内，O_eX_e轴指向格林威治子午线方向，且三轴构成右手坐标系；

3)地理坐标系OX_tY_tZ_t：该坐标系与当地地理经纬度有关，其原点位于载体所处的地球表面，通常选择该坐标系作为导航坐标系，又根据坐标轴定义的不同，可以分为东北天坐标系、北东地坐标系、北西天坐标系等，本发明选取北东地导航坐标系；

4)载体坐标系OX_bY_bZ_b：该坐标系与飞机载体固连，通常取载体重心为原点，OX_b轴沿机身向前，OZ_b轴竖直向上，OY_b轴由右手坐标系确定；

5)星敏感器坐标系OX_sY_sZ_s：该坐标系的原点为星敏感器光学镜头投影中心，OX_s轴、OY_s轴与感光面两边平行，OZ_b轴为沿视轴向外方向，三轴由右手坐标系确定；

6)惯导系统坐标系OX_gY_gZ_g：该坐标系为惯导系统自身坐标系，三轴方向与三个陀螺旋转方向一致，且构成右手坐标系。

(2)相对参考系转换

实际工作中，星敏感器输出姿态数据为星敏感器坐标系OX_sY_sZ_s相对于J2000.0惯性坐标系O_eX_iY_iZ_i的姿态数据，用姿态矩阵Q_i ^s表示；而惯导系统输出姿态数据则为惯导系统坐标系OX_gY_gZ_g相对于地理坐标系OX_tY_tZ_t的姿态数据，用姿态矩阵Q_t ^g表示，显然，星敏感器与惯导系统的相对参考系不同。因此，为了使用高动态星敏感器完成对惯导系统姿态数据的精度校准与评定，首先需要将星敏感器与惯导系统的相对参考系统一，即建立惯性坐标系O_eX_iY_iZ_i与地理坐标系OX_tY_tZ_t的转换关系。如图3所示，可以利用地球坐标系O_eX_eY_eZ_e为过渡坐标系，建立惯性坐标系O_eX_iY_iZ_i与地理坐标系OX_tY_tZ_t的转换关系，有关转换方法在中国专利“一种基于星敏感器确定载体相对于地理坐标系姿态的方法”(CN201110343373.5)中有详细推导。若按照上述专利所提方法求得从惯性坐标系O_eX_iY_iZ_i到地球坐标系O_eX_eY_eZ_e的旋转矩阵为T_i ^e，由地球坐标系到北东地地理坐标系的转换矩阵为T_e ^t，则将星敏感器输出姿态矩阵Q_i ^s转换到地理坐标系下，即星敏感器坐标系相对于地理坐标系的姿态矩阵Q_t ^s如式(1)所示。

$Q_t^s=Q_{ie}^s\·T_e^i\·T_t^e=Q_i^s\·{\left(T_i^e\right)}^T\·{\left(T_e^t\right)}^T---\left(1\right)$

由式(1)可知，变换后星敏感器的姿态数据Q_t ^s与惯导系统的姿态数据Q_t ^g均相对于地理坐标系，即通过上述坐标系变换，可以将星敏感器与惯导系统的相对参考系加以统一。

(3)初始安装误差矩阵求解

当星敏感器与惯导系统的相对参考系统一后，理论上可以利用星敏感器姿态数据对惯导系统数据进行校准与评定。但是实际中，由于安装误差的存在，如图4所示，星敏感器坐标系与惯导坐标系不完全一致，若采用欧拉角表示，则二者的安装误差为(Φ_x,Φ_y,Φ_z)，若采用矩阵表示，则二者的安装误差为B_s ^g。为了便于统一计算，采用矩阵表示方式，即采用B_s ^g表示从星敏感器坐标系旋转至惯导系统坐标系的安装误差矩阵。

由图1可知，星敏感器与惯导系统刚性固连在安装机构上，二者间的误差矩阵B_s ^g为常值偏差矩阵。同时，考虑到惯导系统虽然存在长期漂移，但其短时精度很高，因此静态下、短时间内，利用惯导系统姿态数据与星敏感器姿态数据能够将B_s ^g解算出来，由旋转关系，可以得到B_s ^g的表达式如式(2)所示。

$B_s^g=Q_t^g\·Q_s^t=Q_t^g\·{\left(Q_t^s\right)}^T---\left(2\right)$

实际中，为了提高B_s ^g的估计精度，选取N个采样点，每个采样点采集M对有效数据，则共有N×M个星敏感器与惯导系统的有效数据对，每对数据满足式(3)，式中，表示第i个采样点的第j组对应姿态数据对，B_s ^g为待求的安装误差矩阵。

$Q_{t\left(ij\right)}^g=B_s^g{Q}_{t\left(ij\right)}^s,1\≤i\≤N,1\≤j\≤M---\left(3\right)$

本发明采用QUEST算法求解B_s ^g的最佳估计值，其步骤如下：

①首先，对式(3)进行改写，其等价表达式如式(4)所示，式中，r₁、r₂、r₃和w₁、w₂、w₃分别表示的列向量。

${\left[\begin{array}{ccc}{r}_1& r_2& r_3\end{array}\right]}_{\left(ij\right)}=B_s^g{\left[\begin{array}{ccc}{w}_1& w_2& w_3\end{array}\right]}_{\left(ij\right)},1\≤i\≤N,1\≤j\≤M---\left(4\right)$

进一步的，将式(4)写成列向量分量形式如式(5)所示，式中(r_ijk,w_ijk)表示第i个采样点的第j组对应姿态数据的第k个列向量分量对。

$r_{ijk}=B_s^g{w}_{ijk},1\≤i\≤N,1\≤j\≤M,1\≤k\≤3---\left(5\right)$

②经过上述改写后，式(3)等价为式(5)，而式(5)为典型的Wahba问题，即求解B_s ^g的最佳估计值，等效为求解行列式为1的最优正交矩阵B_s ^g，使得损失函数：

$L\left(B_s^g\right)\≡\frac{1}{2}\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}{\left|r_{ijk}-B_s^g{w}_{ijk}\right|}^2---\left(6\right)$

最小，式中，a_ijk为非负系数，这里取a_ijk恒为1。

③最后，利用QUEST算法，求解B_s ^g的最优估计值。有关QUEST算法的理论推导这里不再展开，仅给出其计算过程及结果如下：

首先，由获得的有效数据，分别计算有关中间变量，其公式如式(7)～(11)所示。

$B=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}{r}_{ijk}{\left(w_{ijk}\right)}^T---\left(7\right)$

$\mathrm{\σ}=tr\left(B\right)=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}{r}_{ijk}\·w_{ijk}---\left(8\right)$

S＝B+B^T(9)

$Z=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}(r_{ijk}\×w_{ijk})---\left(10\right)$

利用上述中间变量，可以将矩阵B_s ^g的最优估计问题，转化为与之等价的四元数的最优估计问题，其公式如(12)所示，式中λ_max为矩阵K的最大特征值，与B_s ^g等价，为λ_max对应的特征向量。

$K{\stackrel{\‾}{q}}_{opt}={\mathrm{\λ}}_{max}{\stackrel{\‾}{q}}_{opt}---\left(12\right)$

对式(12)，在给定精度下可以利用牛顿-拉夫逊迭代法求解则由及式(7)～(10)可得，的最优估计计算过程如式(13)～(16)所示，其中，式(13)中的tr(A)、adj(A)和det(A)分别表示矩阵A的迹、伴随矩阵和行列式。

κ＝tr(adj(S)),△＝det(S)(13)

$\left\{\begin{array}{c}\mathrm{\α}={\left({\widehat{\mathrm{\λ}}}_{max}\right)}^2-{\mathrm{\σ}}^2+\mathrm{\κ}\\ \mathrm{\β}={\widehat{\mathrm{\λ}}}_{max}-\mathrm{\σ}\\ \mathrm{\γ}=({\widehat{\mathrm{\λ}}}_{max}+\mathrm{\σ})\·\mathrm{\α}-\mathrm{\Δ}\end{array}\right.---\left(14\right)$

X＝(αI+βS+S²)Z(15)

${\stackrel{\‾}{q}}_{opt}={\stackrel{\→}{q}}_r+q_0={(q_1,q_2,q_3,q_0)}^T=\frac{1}{\sqrt{{\mathrm{\γ}}^2+{\left|X\right|}^2}}\left\{\begin{array}{c}X\\ \mathrm{\γ}\end{array}\right\}---\left(16\right)$

式(16)所示就是的最优估计，最后，通过四元数与旋转矩阵的转换关系，可以得到B_s ^g的最优估计结果如式(17)所示。

${\hat{B}}_s^g=\left[\begin{array}{ccc}{q_0}^2+{q_1}^2-{q_2}^2-{q_3}^2& 2(q_1{q}_2+q_0{q}_3)& 2(q_1{q}_3-q_0{q}_2)\\ 2(q_1{q}_2-q_0{q}_3)& {q_0}^2-{q_1}^2+{q_2}^2-{q_3}^2& 2(q_3{q}_2+q_0{q}_1)\\ 2(q_1{q}_3+q_0{q}_2)& 2(q_3{q}_2-q_0{q}_1)& {q_0}^2-{q_1}^2-{q_2}^2+{q_3}^2\end{array}\right]---\left(17\right)$

以上就是初始安装误差矩阵B_s ^g的求解过程，需要注意的是，该步骤在整个校准与评定过程中仅需完成一次，即校准与评定前利用上述方法求解后，星敏感器与惯导系统的安装误差保持不变，因此前述求解得到的对于后续校准与评定过程完全适用。

二、数据采集与校准

(1)姿态数据采集

为了实现惯导系统的姿态参数校准，首先需要进行数据采集实验，然后再对实验获得的数据进行校准与评定，而整个惯导系统的姿态数据采集实验过程如图5所示。由图5可知，整个惯导系统的姿态参数校准实验在飞机载体上完成，因此，飞机载体的飞行方案设计至关重要，在飞机具备的机动能力内需要尽可能充分的验证惯导系统的动态性能，同时整个实验应当进行足够长时间，从而考察惯导系统的漂移特性。

(2)姿态数据处理与校准

当惯导系统姿态数据采集实验完成后，需要对实验数据进行校准与评定。如前所述，经过坐标系基准统一后，星敏感器输出姿态数据Q_i ^s能够转换为惯导系统坐标系下的基准姿态数据Q_t ^g _{(s_cal)}，而由式(1)、式(17)可得，该基准姿态数据如式(18)所示，因此，可以利用基准数据Q_t ^g _{(s_cal)}完成对惯导实际输出姿态数据Q_t ^g的精度校准与评定。

$Q_{t(s\_cal)}^g={\hat{B}}_s^g\·Q_i^s\·T_e^i\·T_t^e---\left(18\right)$

实际应用中，姿态矩阵虽然表征了坐标系的旋转关系，具有明确的物理意义，但是姿态矩阵的9个参数(其中3个为独立参数)不具备任何物理意义，因此，姿态矩阵数据不便于直接进行校准与评定。为此，本发明采用与姿态矩阵等价的三轴欧拉角来进行校准与评定，三轴欧拉角为角度量，其表征了坐标系三次旋转的转换关系，具有明确的物理意义，本发明定义欧拉角分别为偏航角、俯仰角和滚转角。

设由星敏感器得到的惯导系统基准姿态欧拉角为该基准值作为惯导系统的姿态真值，若惯导系统测量得到的姿态欧拉角为则惯导系统的绝对测量误差、相对测量误差如式(19)、(20)所示。

但是，由于测量原理不同，现有星敏感器的姿态数据更新率一般低于惯导系统的更新率，且由于时间同步误差，星敏感器与惯导系统的姿态测量时刻不能完全一致，因此，无法直接应用式(19)、(20)对惯导系统姿态测量误差进行计算。为此，需要利用星敏感器所获得的离散姿态真值数据点插值获得姿态真值曲线，然后利用真值曲线对惯导系统姿态参数进行校准与评定。

由于实际飞机载体的姿态连续变化，因此，选用分段三次样条插值方法逼近真值曲线作为较佳实施例，若实验中，星敏感器共采样获得L+1个数据点，则以偏航角数据点(j＝0,1,2,…,L)为例，构造分段三次样条函数：

${\mathrm{\φ}}_s\left(t\right)=\left\{\begin{array}{c}{\mathrm{\φ}}_{s1}\left(t\right),t\∈\[t_0,t_1\]\\ {\mathrm{\φ}}_{s2}\left(t\right),t\∈\[t_1,t_2\]\\ \mathrm{...}\\ {\mathrm{\φ}}_{sL}\left(t\right),t\∈\[t_{L-1},t_L\]\end{array}\right.---\left(21\right)$

满足条件：

①Φ_s(t)在每个小区间[t_j,t_j+1]上为三次多项式；

②Φ_s"(t)在整个区间[t₀,t_L]上连续；

③

④Φ_s"(t₀)＝Φ_s"(t_L)＝0。

则所构造函数Φ_s(t)为偏航角真值曲线，同理构造获得Θ_s(t)、Γ_s(t)分别为俯仰角θ_s、滚转角γ_s的真值曲线，此时应用式(19)、(20)即可对惯导系统姿态测量误差进行计算，同时，还可以对数据标准差进行计算，如式(22)所示，式中，惯导系统共有K+1个测量数据点。

式(19)、(20)和(22)即为惯导系统的误差大小及数据偏差情况的校准与评定公式。

以上就是本发明基于高动态星敏感器的机载惯导系统姿态参数校准方法与装置的技术内容和方案。需要注意的是，以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和范围之内所作的任何修改、等同替换和改进等，均包含在本发明的保护范围之内。

Claims (4)

1.一种机载惯导系统姿态参数校准装置，其特征在于：该校准装置包括高动态星敏感器、GPS同步授时系统、数据采集器、稳压电源和固定支架，其中，固定支架用于固定安装校准装置高动态星敏感器、GPS同步授时系统、数据采集器和稳压电源以及待校准的惯导系统，稳压电源则为高动态星敏感器、GPS同步授时系统和数据采集器以及待校准的惯导系统供电，当需要校准时，将待校准的惯导系统安装在固定支架上，高动态星敏感器将姿态数据及其对应时刻的脉冲信号输出到数据采集器，GPS同步授时系统用于提供标准GPS时间，其中，待校准惯导系统数据直接进行授时，而高动态星敏感器则需要首先对脉冲信号进行授时，然后利用脉冲信号与高动态星敏感器姿态数据一一对应的性质，间接完成对高动态星敏感器输出姿态数据的授时。

2.一种机载惯导系统姿态参数校准方法，利用权利要求1所述的机载惯导系统姿态参数校准装置，其特征在于：该校准方法步骤如下：

步骤一、坐标系基准统一步骤；

将星敏感器与惯导系统的相对参考系统一，即建立惯性坐标系与地理坐标系的转换关系；然后，利用惯导系统姿态数据与星敏感器姿态数据将二者间的常值误差矩阵解算出来；

步骤二、数据采集与校准步骤；

步骤(1)、以恒星为测量对象，通过高动态星敏感器探测天球上的多颗恒星，并进行星点定位、星图识别和姿态解算，提供星敏感器相对于惯性参考系的姿态；以步骤一得到的星敏感器坐标系与惯导系统坐标系的基准统一变换关系，得到惯导坐标系下的飞机姿态基准数据；

步骤(2)、以步骤(1)得到的基准数据为参考，与惯导系统的姿态数据进行比对，从而实现机载环境下对惯导系统姿态精度的动态在线校准。

3.根据权利要求2所述的一种机载惯导系统姿态参数校准方法，其特征在于：步骤一坐标系基准统一的关键在于常值误差矩阵B_s ^g的解算，通过选取N个采样点，每个采样点采集M对有效数据，则共有N×M个星敏感器与惯导系统的有效数据对，每对数据满足式(3)，式中，表示第i个采样点的第j组惯导系统相对于地理坐标系的姿态数据与星敏感器相对于地理坐标系的姿态数据组成的数据对，B_s ^g为待求的常值误差矩阵，

$Q_{t\left(ij\right)}^g=B_s^g{Q}_{t\left(ij\right)}^s,1\≤i\≤N,1\≤j\≤M---\left(3\right)$

然后，采用QUEST算法求解B_s ^g的最佳估计值，其步骤如下：

①首先，对式(3)进行改写，其等价表达式如式(4)所示，式中，r₁、r₂、r₃和w₁、w₂、w₃分别表示的列向量，

${\left[\begin{array}{ccc}{r}_1& r_2& r_3\end{array}\right]}_{\left(ij\right)}=B_s^g{\left[\begin{array}{ccc}{w}_1& w_2& w_3\end{array}\right]}_{\left(ij\right)},1\≤i\≤N,1\≤j\≤M---\left(4\right)$

进一步的，将式(4)写成列向量分量形式如式(5)所示，式中(r_ijk,w_ijk)表示第i个采样点的第j组对应姿态数据的第k个列向量分量对，

$r_{ijk}=B_s^g{w}_{ijk},1\≤i\≤N,1\≤j\≤M,1\≤k\≤3---\left(5\right)$

②经过上述改写后，式(3)等价为式(5)，而式(5)为典型的Wahba问题，即求解B_s ^g的最佳估计值，等效为求解行列式为1的最优正交矩阵B_s ^g，使得损失函数：

$L\left(B_s^g\right)\≡\frac{1}{2}\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}|r_{ijk}-B_s^g{w}_{ijk}{|}^2---\left(6\right)$

最小，式中，a_ijk为非负系数，这里取a_ijk恒为1；

③最后，利用QUEST算法，求解B_s ^g的最优估计值：

首先，由获得的有效数据，分别计算有关中间变量，其公式如式(7)～(11)所示：

$B=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}{r}_{ijk}{\left(w_{ijk}\right)}^T---\left(7\right)$

$\mathrm{\σ}=tr\left(B\right)=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}{r}_{ijk}\·w_{ijk}---\left(8\right)$

S＝B+B^T(9)

$Z=\underset{i=1}{\overset{N}{\Σ}}\underset{j=1}{\overset{M}{\Σ}}\underset{k=1}{\overset{3}{\Σ}}{a}_{ijk}(r_{ijk}\×w_{ijk})---\left(10\right)$

利用上述中间变量，可以将矩阵B_s ^g的最优估计问题，转化为与之等价的四元数的最优估计问题，其公式如(12)所示，式中λ_max为矩阵K的最大特征值，与B_s ^g等价，为λ_max对应的特征向量，

$K{\stackrel{\‾}{q}}_{opt}={\mathrm{\λ}}_{max}{\stackrel{\‾}{q}}_{opt}---\left(12\right)$

对式(12)，在给定精度下可以利用牛顿-拉夫逊迭代法求解则由及式(7)～(10)可得，的最优估计计算过程如式(13)～(16)所示，其中，式(13)中的tr(A)、adj(A)和det(A)分别表示矩阵A的迹、伴随矩阵和行列式，

κ＝tr(adj(S)),Δ＝det(S)(13)

$\left\{\begin{array}{c}\mathrm{\α}={\left({\widehat{\mathrm{\λ}}}_{max}\right)}^2-{\mathrm{\σ}}^2+\mathrm{\κ}\\ \mathrm{\β}={\widehat{\mathrm{\λ}}}_{max}-\mathrm{\σ}\\ \mathrm{\γ}=({\widehat{\mathrm{\λ}}}_{max}+\mathrm{\σ})\·\mathrm{\α}-\mathrm{\Δ}\end{array}\right.---\left(14\right)$

X＝(αI+βS+S²)Z(15)

${\stackrel{\‾}{q}}_{opt}={\stackrel{\→}{q}}_r+q_0={(q_1,q_2,q_3,q_0)}^T=\frac{1}{\sqrt{{\mathrm{\γ}}^2+|X{|}^2}}\left\{\begin{array}{c}X\\ \mathrm{\γ}\end{array}\right\}---\left(16\right)$

式(16)所示就是的最优估计，最后，通过四元数与旋转矩阵的转换关系，可以得到B_s ^g的最优估计结果如式(17)所示，

${\hat{B}}_s^g=\left[\begin{array}{ccc}{q_0}^2+{q_1}^2-{q_2}^2-{q_3}^2& 2\left(q_1{q}_2+q_0{q}_3\right)& 2\left(q_1{q}_3-q_0{q}_2\right)\\ 2\left(q_1{q}_2-q_0{q}_3\right)& {q_0}^2-{q_1}^2+{q_2}^2-{q_3}^2& 2\left(q_3{q}_2+q_0{q}_1\right)\\ 2\left(q_1{q}_3+q_0{q}_2\right)& 2\left(q_3{q}_2-q_0{q}_1\right)& {q_0}^2-{q_1}^2-{q_2}^2+{q_3}^2\end{array}\right]---\left(17\right).$

4.根据权利要求2所述的一种机载惯导系统姿态参数校准方法，其特征在于：步骤二数据采集与校准，由于测量原理不同，现有星敏感器的姿态数据更新率一般低于惯导系统的更新率，且由于时间同步误差，星敏感器与惯导系统的姿态测量时刻不能完全一致，故而无法直接应用星敏感器数据对惯导系统姿态参数进行校准，而是需要利用星敏感器所获得的离散姿态真值数据点插值获得姿态真值曲线，然后利用真值曲线对惯导系统姿态参数进行校准与评定。