基于误差时间特性的INS/GNSS观测噪声方差在线计算方法

邵梦晗,高晓颖,吕建强

(1.北京航天自动控制研究所， 北京 100854； 2.宇航智能控制技术国家级重点实验室， 北京 100854)

精确制导与精确入轨效能主要受制于导航系统性能，导航系统的精度从根本上决定了后续制导、控制的精度。目前工程上广泛使用基于卡尔曼滤波的INS/GNSS组合导航，其滤波效果依赖于噪声统计特性等先验信息。但对飞行环境复杂、高机动的载体，多路径效应、可见卫星数目及几何分布、外界电磁波干扰等多种因素并存，造成GNSS观测噪声多变且难以预测，为此涌现了多种自适应滤波算法[1-4]。但大部分自适应算法都是从卡尔曼滤波算法本身出发，获得观测量与滤波参数的自适应调节关系，同时又引入了其他需要事先确定的参数，导致其对不同噪声环境的适应性、可靠性不高，难以在实际工程中应用。文献[5]提出一种组合导航滤波参数在线估计的并行架构，并将遗传算法引入进行参数优化，但是适值函数难以确定。文献[6]提出了一种根据传感器不同测量特性的自适应滤波算法，简单易行，但其只进行了理论分析，且仅限于组合导航开环校正使用。然而目前工程上广泛使用速度、位置反馈校正，甚至全状态量闭环校正方式。为此，本文针对组合导航闭环校正方式，依据常用惯性导航短时高精度和GNSS测量噪声的时间不相关特性，提出一种基于小样本统计的在线计算GNSS观测噪声方差阵的方法；并利用实测数据验证了这种误差特性，为该思路的可行性提供了数据支撑；最后给出了具体的实现步骤，并利用一组飞行仿真数据对该方法进行了验证。

1 闭环校正下GNSS观测噪声方差在线计算方法

1.1 发射惯性系下组合滤波模型

在发射惯性系下，以速度、位置、姿态、加速度计零偏、陀螺零偏的误差量作为组合滤波状态量X(t)，以INS与GNSS的位置差、速度差作为观测量Z(t)，构建如下滤波模型：

(1)

Z(t)=H(t)X(t)+Wgnss(t)

(2)

X(t)=(δvx,δvy,δvz,δrx,δry,δrz,

δφx,δφy,δφz,δK0x,δK0y,δK0z,δD0x,δD0y,δD0z)T

Z(t)=(δvx,δvy,δvz,δrx,δry,δrz)T

其中，δvx,δvy,δvz为速度误差；δrx,δry,δrz为位置误差；δφx,δφy,δφz为姿态角误差；δK0x,δK0y,δK0z为加速度零次项误差；δD0x,δD0y,δD0z为陀螺零次项误差；Wins(t)、Wgnss(t)分别是系统白噪声阵、观测白噪声阵；F(t)是系统矩阵；G(t)是系统噪声驱动阵；H(t)是观测矩阵。

1.2 惯性导航短期高精度信息利用

卫星导航速度误差在0.1 m/s量级、位置误差在10 m量级，且测量误差与时间不相关，误差不随着时间积累；常用激光惯组陀螺零偏在0.01°/h量级、加速度计零偏在10-5g量级，惯性导航更新频率快，短期精度高，但误差随着时间累积而发散。但是惯性导航在1 s、5 s乃至更长一段时间内的测量误差都是小于卫星导航测量误差的，为此，考虑充分挖掘惯性导航的短期高精度信息，对GNSS观测噪声进行在线估计。但载体实际飞行时无真实值参考，因此，联想到利用INS、GNSS各自相邻时刻增量信息：惯组测量噪声远小于GNSS接收机测量噪声，且相邻时刻作差能够抵消一部分惯组不随时间变化的误差，导致纯惯性导航短期增量精度高于卫星导航的增量精度，以此来估计GNSS观测噪声方差阵。

1.3 GNSS观测噪声方差在线计算方法

以误差量为卡尔曼滤波器状态量，得到的估计有两种利用方法[7-8]：一种是将估计作为惯导系统输出的校正量，称为开环法；另一种是将估计反馈到惯导系统中，作为惯性导航力学编排方程中的相应参数或校正量，即闭环法。开环校正的滤波方程就是离散化后的卡尔曼滤波方程；闭环校正每次校正后需要把状态量清零，其余与开环校正无异。

开环校正下，INS、GNSS的k时刻速度误差、位置误差状态估计量：

(3)

其中，Xv,r(k)是速度误差、位置误差状态量真值；fins(k)是惯性导航误差趋势项。

闭环校正下，组合滤波反馈修正公式：

(4)

由惯性导航更新公式：

(5)

(6)

其中，dδXv,r(k-1)是速度、位置修正量误差。取GNSS更新周期为组合滤波周期，分别求INS、GNSS相邻滤波时刻的状态增量：

ΔINS(k)=Xins(k)-Xins(k-1)=[X(k)-X(k-1)]+

[fins(k-1)-fins(k)]+[Wins(k-1)-Wins(k)]+

[dδXv,r(k-2)-dδXv,r(k-1)]

(7)

ΔGNSS(k)=Xgnss(k)-Xgnss(k-1)=

[X(k)-X(k-1)]+[Wgnss(k-1)-Wgnss(k)]

(8)

式(7)、式(8)作差：

ΔINS(k)-ΔGNSS(k)=[fins(k-1)-fins(k)]+

[Wins(k-1)-Wins(k)]-[Wgnss(k-1)-Wgnss(k)]+

[dδXv,r(k-2)-dδXv,r(k-1)]

(9)

求式(9)的自协方差，并考虑Wins(k)、Wgnss(k)是零均值白噪声，有：

[fINS(k-1)-fINS(k)]2+Q(k-1)+Q(k)+R(k-1)+

R(k)+Pv,r(k-2)+Pv,r(k-1)

(10)

其中，Q(k)是系统噪声方差阵，R(k)是观测噪声方差阵，Pv,r(k-2)、Pv,r(k-1)分别是k-2、k-1时刻组合滤波速度、位置状态量误差方差阵，是6×6维矩阵。

考虑惯组测量噪声远小于GNSS接收机测量噪声Q(k)<

[fINS(k-1)-fINS(k)]2<

(11)

故上式(10)可写为

R(k-1)+R(k)+Pv,r(k-2)+Pv,r(k-1)

(12)

从而，可得到GNSS观测噪声方差阵R(k)的统计估计：

[Pv,r(k-2)+Pv,r(k-1)]}

(13)

2 实测数据验证与定量分析

截取某飞行器某段时间的纯惯性导航、卫星导航的实测数据，激光惯组和GNSS接收机均独立、正常工作。惯性导航更新周期是0.02 s，卫星导航更新周期是1 s。为同步求增量，惯性导航需要进行1 s的递推解算，然后分别求INS、GNSS相邻1 s的速度增量Δv、位置增量Δr，绘制出惯性导航与卫星导航的增量曲线如图1、图2所示。

由图1、图2可以看出：GNSS结果增量在惯性导航结果增量附近波动。两者的增量均值误差可忽略不计，任取某时间段的样本数据，计算纯惯性导航速度增量标准差为0.009 6 m/s，位置增量标准差为0.566 1 m；GNSS速度增量标准差为0.198 4 m/s，位置增量标准差为6.236 7 m。可见，纯惯性导航增量的标准差远小于对应的GNSS结果增量的标准差，即：

[fINS(k-1)-fINS(k)]2+Q(k)+Q(k-1)<

(14)

进而说明式(12)取近似是合理的。

图1 速度增量曲线

图2 位置增量曲线

进一步，根据实际飞行情况进行定量分析，假设载体运动速度v=3 000 m/s，加速度a=20 m/s2，速度方向误差增量Δθ=1′，则Δt=1 s时的速度误差增量、位置误差增量分别为

Δvins=a×Δθ×Δt=20×1/60/180×π=0.005 8 m/s

Δrins=v×Δθ×Δt=3 000×1×1/60/180×π=

0.872 639 m

可以看出：即使在较高速度和大误差角增量的假设下，1 s内纯惯性导航的速度误差增量只有0.006 m/s、位置误差增量只有0.9 m。而GNSS测速误差一般为0.1 m/s量级，定位误差一般为10 m量级。因此，1 s内惯性导航误差增量的平方值远小于GNSS定位误差的方差[6]，所以式(11)的假设是合理的，因此式(12)是符合客观实际的。

3 算法实现

考虑实际算法的稳定性和可行性，具体实现步骤如下：

1) 计算INS、GNSS各自相邻时刻速度、位置状态增量ΔINS、ΔGNSS，并计算：

B(k)=ΔINS(k)-ΔGNSS(k)

(15)

B(k)序列含有GNSS测量噪声信息，因此可以通过对一段时间的样本数据进行统计，得到观测噪声阵的估计。样本容量越大，统计结果的可靠性越高，但当实际噪声发生变化时，样本容量越大，这种动态信息被平滑，即当前时刻的动态信息淹没在大量的历史信息中，导致统计结果不能及时跟踪动态的噪声信息，实时性较差。样本容量越小，统计结果的可靠性降低。所以，通过滑动窗口利用最新一段导航信息进行方差统计，具体公式为

[Pv,r(k-2)+Pv,r(k-1)]}

(16)

其中，k为当前滤波时刻，k>M；M为滑动窗口长度，一般为20～50 s，可根据实际飞行情况选取。

2) 粗差剔除。B(k)受传感器观测粗差的影响较大，因此在计算前需要剔除粗差。一般认为，高精度激光惯组出现粗差的可能性很小，因此，主要是考虑GNSS出现粗差的情况。本文采用标准化新息作为假设检验统计量进行粗差探测，具体方法可参考文献[9-11]。若探测到某时刻存在观测粗差，将该时刻数据剔除，不参与统计估计。

3) 观测噪声方差阵R(k)的估计。为避免出现小样本统计下R(k)变化较为剧烈的情况，引入指数加权法[8]更新R(k)阵：

(17)

其中，b为遗忘因子，一般取0.95～0.995。

4 仿真验证

4.1 仿真设置

采用一组模拟飞行数据进行仿真验证，在发射惯性系下，初始速度为[406,0.1,-53]m/s；初始位置为[0.1,0.1,0.1]m；初始姿态为[90,-0.1,0.1]°；陀螺零偏0.1°/h，测量噪声0.05°/h；加速度计零偏0.002g，测量噪声50 μg；组合滤波时间为1 302 s。仿真设置R(k)阵为对角阵，其值是GNSS观测白噪声的方差。为考察闭环校正下在线计算R(k)阵的效果，在280～580 s设置观测噪声方差发生突变，具体R(k)阵对角线元素仿真真值如表1所示。

表1 R(k)仿真真值

为进一步验证组合导航闭环校正下在线计算方法的有效性，将其与同条件下标准卡尔曼滤波、理想方法对比分析。其中，标准卡尔曼滤波方法与在线计算方法的观测噪声方差阵初值相同，为充分验证该方法的有效性，设置观测噪声方差阵初值为：Rv0=0.04(m/s)2，Rr0=400(m)2，与仿真真值相差较大；理想方法下观测噪声方差阵始终等于仿真真值，其结果可以作为参考值对比。

4.2 仿真结果与分析

将在线计算方法、标准卡尔曼滤波方法、理想方法解算的速度、位置分别与仿真轨迹真值作差，得到速度误差、位置误差如图3、图4所示。

图3 速度误差

图4 位置误差

同时统计了观测噪声发生突变时间段的速度合成误差、位置合成误差均值和均方根，如表2所示。

分析图3、图4，在初始观测噪声方差阵误差较大情况下，本文提出的在线计算方法能够较好地根据实际情况进行调整，并能及时跟踪噪声突变，精度与理想方法较为接近；而标准卡尔曼滤波算法始终采用较差的初始观测噪声方差阵计算，速度、位置估计出现了较大的误差。分析表2，噪声突变段，在线计算方法解算的速度误差均值0.016 m/s，位置误差均值0.567 m，速度误差均方根0.097 m/s，位置误差均方根1.255 m，精度与理想方法相当；而标准卡尔曼滤波得到的速度误差均值0.1 m/s，位置误差均值13.864 m，速度误差均方根0.210 m/s，位置误差均方根9.416 m。可见，本文提出的在线计算滤波方法性能明显优于标准卡尔曼滤波算法。

表2 噪声突变时间段(280～580 s)导航误差均值、均方根

注：表2中方法1是本文提出的在线计算观测噪声方差阵方法；方法2是标准卡尔曼滤波方法；方法3是理想方法。

绘制3种方法下的每滤波时刻R(k)阵对角线元素值如图5所示。

图5 速度、位置观测噪声方差

可以看出：标准卡尔曼滤波方法测量噪声方差始终为初始值，误差较大；而在线计算R(k)方法是以过去一段时间的滑动窗口数据为样本进行统计计算，同时为避免大幅震荡采用了指数加权，因此存在一定的振荡和滞后，但是能够随着真实观测噪声的变化而变化，始终跟踪实际值，效果较好。

5 结论

本文提出了一种基于INS、GNSS不同误差时间特性的在线计算GNSS观测噪声方差阵的方法。该方法利用滑动时间窗口内的增量样本数据在线计算观测噪声方差阵，避免了传统自适应滤波算法引入其他参数调整的问题，对工程上广泛使用的组合导航闭环校正适用。仿真试验表明该方法能够准确跟踪观测噪声的变化，获得较高精度的导航状态量，对环境适应性较强，可应用于飞行环境复杂、高机动载体的导航测量中。