基于最大熵卡尔曼滤波的机载SINS星敏感器辅助空中标定技术

王 瑞，张志亮，郑百东，戴洪德

(海军航空大学 航空基础学院，山东 烟台 264001)

0 引言

由于惯性器件性能参数及系统特性随着时间变化，惯性导航系统的精度会随时间降低，在一定时间之后，需要对惯性导航系统的误差进行标定，以便进行补偿从而提高导航精度[1]。惯性导航系统标定方法根据标定时是直接采用惯性器件输出参数还是导航参数可以分为分立式标定[2]和系统级标定[3]。分立式标定需要将惯性导航系统从飞机等载体上卸载，并通过转台精准模拟飞机运动，对惯性导航系统误差进行标定补偿。分立式标定的精度较高，但是需要繁琐的装卸操纵、特殊的转台设备，且时间效率差[4]，影响正常使用。系统级标定则可以直接在载体上在线完成，通过对比惯性导航系统输出的导航参数和其他系统较为准确的导航参数，得出差值作为观测量，通过相应的滤波算法估计惯性导航系统的误差，对惯性导航系统进行补偿与标定。系统级标定具有方便、快捷、效率高的特点，已成为惯性导航系统标定方法的主要研究热点。

在系统级标定中，根据借助的外部导航参数不同，惯性导航系统能够标定的参数的效果是不同的。利用里程计可采用位置误差进行标定[5-6]，利用GPS信号可采用速度误差进行标定[7]，利用星敏感器可采用姿态误差进行标定[8]。由于位置是由速度积分而来，且位置信息的可观测性较差，利用速度误差和姿态误差对惯性导航系统进行标定的效果优于位置误差的标定效果。单独采用一项外部设备标定惯性导航系统有可能造成其他误差不可观测，待标度误差的估计收敛速度也较慢[9]。在滤波算法中，卡尔曼滤波(KF，Kalman filter)对误差的估计效果较好，但非高斯噪声的存在限制了卡尔曼滤波对误差估计精度的提升。

本文针对机载捷联惯性导航系统(SINS, strapdown inertial navigation system)空中标定问题，结合星敏感器和GPS差分信号，对比得到姿态偏差和速度偏差，以其作为观测量，并用最大熵卡尔曼滤波(MCKF, maximum correntropy Kalman filter)方法解决观测值中存在非高斯噪声的问题，以提升惯性导航系统误差标定的速度和精度。

1 星敏感器辅助SINS自标定技术

1.1 总体方案

星敏感器是天文定姿设备，通过拍摄天空恒星图像，与存储的星图对比，确定飞机姿态，其精度较高，且误差不累积。能够输出高精度的姿态信息，且不随时间漂移，自主性强。将星敏感器输出的姿态信息与SINS的姿态信息对比，得出姿态偏差。为提高系统整体可观测性，同时引入GPS差分信号，与SINS输出的速度信号对比，得出速度偏差。以姿态偏差和速度偏差作为系统级标定模型的观测值，经由卡尔曼滤波器估计出系统真实误差项，对SINS进行标定，整体方案如图1所示。为增加姿态可观测度，设定飞机进行机动飞行。

图1 星敏感器辅助惯性导航系统自标定方框图

1.2 星敏感器工作原理

星敏感器主要由3部分构成，分别为成像器件、图像传感器和数据解算电路，涉及到光学、机械、电子等领域。按照其工作原理进行划分，又可分为成像系统和图像处理系统。星敏感器工作流程如图2所示，首先将光学镜头朝向天空，对镜头范围内的恒星进行成像处理，成像于图像敏感器芯片上，将光学信息转换为电信号，经支持电路将电信号放大为电压模拟信号，随后由A/D转换器转换为数字信号，此时，星敏感器采集的光学成像成功转换为可以存储的数字信号，且存在一一对应的关系。然后，数据处理系统扫描上面存储的数字星图，找出其中的恒星，并计算出观测到恒星的亮度信息和在像平面上的位置信息。最后，经星图识别确定星图内恒星与导航星库中恒星的对应关系，通过姿态解算方法确定星敏感器的三轴姿态，进而为载体提供高精度姿态信息，实现其导航的目的。

图2 星敏感器工作流程框图

其中，星敏感器成像原理如图3所示，O-xsyszs代表星敏感器坐标系，O-xvyvzv代表CCD面阵坐标系，光轴OOs与zv轴、zs轴均重合，f代表光学透镜的焦距。p点为恒星于CCD面阵上的点，其坐标为(x,y)，此时，星光矢量在O-xsyszs中的方向向量为：

(1)

式中，Vs代表观测噪声。

图3 星敏感器成像原理图

假设星敏感器镜头内能够观测到n颗恒星，根据上式可得星光矢量在O-xsyszs中的方向向量分别为(xs1,ys1,zs1)(xs2,ys2,zs2)(xs3,ys3,zs3)…(xsn,ysn,zsn)，搜索导航星库后可知，这n颗恒星在惯性空间中方位分别为(xi1,yi1,zi1)(xi2,yi2,zi2)(xi3,yi3,zi3)…(xin,yin,zin)，则星敏感器姿态转移矩阵满足下式：

(2)

也可写为:

S=GA

(3)

式(2)和式(3)相互对应，其中A为载体的姿态转移矩阵。

当n=3时，

A=G-1S

(4)

当n>3时，根据最小二乘法原理得：

A=(GTG)-1(GTS)

(5)

由上述分析可知，在惯性空间内，恒星相对于地球的方位始终保持不变，星敏感器正是以此来计算载体的高精度姿态信息。星敏感器的测量精度会受到地球的章动、自身误差等影响，但总体影响较小，远小于1″，也不会随时间不断积累。

综上，星敏感器可以为本文自标定系统提供姿态信息，从而提高机载SINS自标定系统的精度。

1.3 陀螺与加速度计误差模型

机载SINS自标定为现场在线标定，没有拆卸操作，不影响安装误差。在代表机体坐标系b和惯性坐标系i下，陀螺仪误差模型[10]可以用式(6)表示:

(6)

同样在机体坐标系下，加速度计误差模型如式(7)所示:

δfb=δKafb+▽a

(7)

1.4 SINS误差方程

取n为导航坐标系，机载接连惯性导航系统的姿态误差方程如式(8)所示，速度误差方程如式(9)所示，位置误差方程如式(10)～(12)所示:

(8)

(9)

(10)

(11)

(12)

1.5 状态方程

根据SINS误差方程，可以构建系统自标定误差模型，但由于不进行惯导及惯导内部器件的重新拆装，即可认为安装误差不发生改变。所以主要对陀螺漂移和加速度计零偏，及相应的标度因数误差进行标定。捷联惯性导航系统在标定时，需对姿态角误差φE、φN、φU，速度误差δVE、δVN，位置误差δL、δλ、δh，陀螺仪漂移、陀螺仪标度因数误差、加速度计零偏、加速度计标度因数误差等进行估计与补偿，将以上误差项作为自标定模型的状态变量，构建系统状态方程如式(13)所示:

(13)

式中，状态变量为X(t)。

X(t)=[φEφNφUδVEδVNδVUδLδλδhεbx

εbyεbz▽ax▽ay▽azδKaxδKayδKazδKgxδKgy

δKgz]T，系统随机噪声W(t)=[wgxwgywgzwaxwaywaz]T，由陀螺仪随机噪声wgxwgywgz和加速度计随机噪声waxwaywaz两部分组成。F(t)和G(t)分别为系统矩阵和噪声输入矩阵:

(14)

式中，FN为9×9维惯导系统误差矩阵，012×21为元素为0的12×21阶矩阵。FS和FM分别如下所示：

其中:元素矩阵F1和F2是陀螺仪和加速度计输出在导航坐标系的投影。

式中，C11、C12、C13等是载体系到导航系转移矩阵的元素。

1.6 量测方程

采用GPS辅助SINS空中自标定方法仅以速度误差作为量测量，标定天向陀螺时漂移估计较慢，效果较差。星敏感器能够求出高精度姿态信息，但对加速度计误差不可观测。所以，在此处将星敏感器与SINS的姿态偏差和差分GPS与SINS的速度偏差，作为自标定模型的观测值量，构成组合量测模型，其量测方程如式(15)所示:

Z(t)=HX(t)+V(t)

(15)

式中，观测量Zav(t)=[φEφNφUδVEδVN]T，Hav=[I5×505×16]5×21，V(t)为量测噪声，其中I5×5为5×5阶单位阵。可以看出，量测方程采用的是“姿态+速度”匹配量测模型，与位置匹配两侧模型相比，可以提升自标定的速度。

2 基于最大熵的卡尔曼滤波

2.1 最大熵准则

相关熵是一种两个随机变量间的广义相似度量，对于联合分布函数为FXY(x,y)的两个随机变量X和Y是实域向量，其相关熵定义如下[12]：

(16)

式中，E[κ(X,Y)]表示κ(X,Y)的期望，κ(.,.)表示Mercer核函数，此处取高斯核函数，即:

(17)

式中，e=x-y且σ>0表示核函数的带宽[13]。

一般情况下，由于已知数据有限且FXY(x,y)未知，可由样本平均估计值计算相关熵：

(18)

由估算值可知，V(X,Y)为有界正数，当X=Y时，V(X,Y)取得最大值。对高斯κ(.,.)泰勒级数展开[14]：

(19)

由式(19)可以看出，变量X和Y的相关熵是对X-Y的偶次幂在取均值后，再加权求和。相关熵中有高阶矩信息，若噪声为非高斯分布，最大相关熵的鲁棒性可显著改善。

2.2 最大熵卡尔曼滤波

基于最小均方差准则的经典卡尔曼滤波器只能精确到误差的二阶统计特性，但是如采用最大相关熵准则作为卡尔曼滤波器的滤波准则，则可以应用误差的高阶统计特性，从而改善非高斯噪声条件下的滤波器鲁棒性。

SINS标定模型离散化后，式(18)和式(20)联立改写为：

(15)

(21)

式中，I为单位阵，

(22)

Dk=WkXk+Ek

(23)

由于Ek为白噪声，取损失函数JL为：

(24)

式中，di,k是Dk的第i行元素，L=n+m是Dk的维数。

由最大相关熵准则可知[16]，Xk的最优估计为使损失函数JL达到最大时的值：

(25)

式中，ei,k=di,k-wi,kXk。

损失函数的最大值可通过求偏导数为0的Xk：

(26)

化简得：

(27)

代入ei,k得：

(28)

以矩阵形式表示：

(29)

其中:

(30)

Cx,k= diag(Gσe1,k,…,Gσen,k)

(31)

Cy,k= diag} (Gσen + 1,k,…,Gσen + m,k)

(32)

将式(29)化简后得[14]：

(33)

式中，

(34)

由此整理可得最大熵卡尔曼滤波器的算法流程：

4)更新增益矩阵;

5)状态更新;

6)协方差阵更新;

由此，最大相关熵卡尔曼滤波通过高阶矩信息，提高了非高斯噪声条件下的鲁棒性，可以解决机载捷联惯导系统在空中标定时难以解决的非高斯噪声问题。

3 仿真实验与分析

3.1 仿真条件

图4 SINS解算原理图

SINS输入的加速度计和陀螺仪数据与飞机运动模型相一致，并叠加各误差项，为自标定算法提供数据支撑。

3.2 仿真方法

不同于传统标定方法，通过采用转台转动激励出SINS相关误差参数，实现惯导误差参数的标定，本文基于常规起飞过程的机载SINS自标定方法，仅利用飞机常规起飞过程中的机动，激励出惯性导航系统内部惯性器件的误差，依据系统级标定原理，以机场跑道附近的高精度差分GPS和惯性导航输出速度差值作为观测量，应用Kalman滤波对惯导系统的状态量进行估计，实现机载SINS的误差标定。

通过飞行轨迹的设计，能够使得系统更好地激励出系统状态参数，从而使机载SINS元器件误差能被更好的标定出来。

由于在飞行过程中GPS存在精度要求无法长时间保证的问题，则充分利用起飞阶段飞机的机动以及机场跑道附件能够布设高精度GPS差分站的优势，将空中自标定的轨迹设置在机场附近。结合常规起飞阶段的仿真分析，设计恰当的飞行轨迹，且保证充分利用机场附近高精度的差分GPS。根据某型飞机的起飞机动特性以及机场的数据，在尽可能不增加飞行员额外操作的原则下设计起飞标定动作共5个阶段，将起飞阶段和盘旋阶段相结合，共588 s。

飞行轨迹如图5所示，假设飞机位于5 km的高度，初始速度为200 m/s，常规的机动动作主要分为5个阶段：

1)匀速阶段：飞机保持初始状态飞行。

2)加速阶段：飞机以2.5 m/s2加速度将飞机速度从200 m/s提升到275 m/s。

3)爬升阶段：飞机在机场上方5 km高度处，以275 m/s速度进行爬升，将飞机以2/s的俯仰角拉升至30°；随后保持角度爬升一段时间；最后将飞机改平，爬升阶段完成。

4)滚转阶段：飞机在机场上方5 km高度处，以275 m/s速度进行滚转，首先以5/s左滚转到45，随后保持该状态飞行一段时间，最后以5/s右滚转将飞机调整到水平，保持原状态继续飞行。

5)转弯阶段：飞机在机场上方5 km高度处，以275 m/s速度进行转弯，首先以9/s右转弯到90，随后保持该状态飞行一段时间，最后以9/s左转弯将飞机调整到原方向，保持原状态继续飞行。

通过飞机机动激励SINS输出参数信息，利用最大相关熵卡尔曼滤波器对各项误差参数进行标定，验证设计的方法有效性。

3.3 仿真参数

机载捷联惯导系统的仿真参数如表1所示。

表1 仿真参数设定

星敏感器的姿态角误差φE=10″，φN=10″，φU=20′，GPS位置误差1 m，速度误差0.2 m/s。

基于以上误差参数，设定协方差初始阵为：

P=diag([10″；10″；20′；0.2 m/s；0.2 m/s；0.2 m/s；1 m；1 m；1 m；0.03°/h；0.03°/h；0.03°/h；100 μg；100 μg；100 μg；400×10-6；400×10-6；400×10-6；400×10-6；400×10-6；400×10-6]×10)2

噪声方差初始矩阵为：

由于采用了星敏感器和GPS信号作为空中标定的辅助信号，所以在观测量中引入姿态和航向三轴角度信号和水平两个速度信号。在非高斯噪声条件下，设定核函数的带宽σ=0.8，量测噪声为:

仿真时，飞机航迹由爬升、转弯、直线飞行等段组成，如图5所示。飞机的初始状态为：高度5 km，速度200 m/s，保持匀速直线飞行，以黑色原点代表初始位置。之后，进入分别进行爬升、转弯、滚转等机动动作，仿真时间共588 s。分别采用最大相关熵卡尔曼滤波和经典卡尔曼滤波依据GPS信号和星敏感器信号，对机载捷联惯导系统进行标定。仿真结果如图6和图7所示，在图中，虚线表示采用经典卡尔曼滤波后惯性器件的误差估计曲线，实线表示采用最大熵卡尔曼滤波后惯性器件的误差估计曲线，直实线代表各误差项的设定值。

图5 飞机仿真航迹

图6 陀螺仪漂移ε与加速度计零偏▽估计曲线

图7 加速度计标度因数误差δKa与陀螺仪标度因数δKg误差估计曲线

3.4 仿真结果与分析

由图6可以看出，采用了最大熵卡尔曼滤波的星敏感器辅助SINS自标定，在三轴陀螺漂移的估计中均优于经典卡尔曼滤波，MCKF的收敛速度和估计精度均高于经典卡尔曼滤波。尤其是在z轴陀螺仪漂移的估计上，MCKF的精度显著高于经典卡尔曼滤波，并且克服了无星敏感器辅助时，z轴陀螺仪漂移估计发散的问题[17]。但在加速度计零偏的估计上，MCKF的优势较小，尤其是在z轴上的估计，难以克服惯导系统高度通道发散的问题。但是，对于SINS垂直通道，通常方法是借助于大气数据参数来修正高度，因此，在后续研究中，应借助大气数据系统的高度参数，提高对垂直通道的观测度和实现标定效果的提升。同时对比文献[17]可以看出，星敏感器对加速度计零偏估计的影响较小，这是由于星敏感器主要用于测量姿态数据，其对加速度的检测影响较小。

由图7可以看出，MCKF算法对加速度计标度因数误差的估计精度比较高，但是收敛速度与KF相当；对陀螺仪标度因数误差的估计精度和收敛速度优于KF。同时，和不采用星敏感器的自标定技术相比[17]，估计精度略有提升。标度因数误差反应惯性测量器件的一致性，反应的是器件本身的测量特性，非高斯噪声对标定因数误差影响较小。因此，可以看出，不管是否采用最大熵原则，对标定因数的标定的影响较小。

不同滤波技术下，陀螺仪漂移、加速度计零偏以及相应的标度因数误差的标定结果如表2所示。

表2 非高斯条件下不同滤波技术的误差标定结果

由表2可知，在非高斯噪声条件下，MCKF的标定结果明显优于经典卡尔曼滤波。但是对于MCKF条件下，星敏感器的采用对标度因数的误差标定提升不够明显。

4 结束语

对于机载SINS的空中自标定，为提升陀螺仪漂移和加速度计零偏的标定精度，结合前期的SINS误差模型的可观性分析，分别采用了星敏感器和GPS信号辅助机载SINS的空中自标定。由于空中飞行时，飞机机动、机体振动等因素引起的附加噪声是非高斯噪声，在此环境下，为消除非高斯噪声对自标定精度的影响，采用了最大相关熵卡尔曼滤波方法。在最大相关熵原则中，引入了高阶矩信息，可以提升滤波器的鲁棒性，从而显著提升了机载SINS误差的自标定精度和速度。为验证所提自标定方法的有效性，以SINS模型和飞机运动模型为基础，通过设计飞机机动航迹，激励SINS各项误差，采用最大相关熵卡尔曼滤波算法来估计误差项。同时和未采用星敏感器辅助的SINS空中标定得到的误差项进行了对比。通过对比分析可以看出，在星敏感器辅助下，最大相关熵卡尔曼滤波算法对陀螺仪漂移和加速度零偏的标定精度有了显著提升，尤其是天向陀螺仪的漂移，标定速度和精度都得到了提升。但是高度通道存在发散情况，因此需借助其他方法提升垂直加速度计零偏的标定精度。同时，本文在建立误差模型时，仅对惯性元器件的常值误差、刻度因数误差和安装误差进行分析讨论，未对惯性元器件温度、高阶项误差等进行深入分析，在后续工作中可建立更全面的误差模型，进一步提高误差项的标定精度。