三段判别域与最小二乘拟合的抗差滤波算法

蔡保杰, 邵 雷

(空军工程大学防空反导学院, 陕西 西安 710051)

0 引 言

导航系统与日常生活关系越来越密切,无论是车载导航还是飞机、航天器等都需要精确的定位。当行驶到隧道或高楼林立的街道,或是在现代战争中卫星信号被大功率杂波干扰,都会使卫星导航观测值出现故障,对导航精度造成严重影响[1-2],抗差滤波便是在这样的背景下被提出的,Krarup[3]等人提出理论,Caspary[4]完善了这一理论。

抗差滤波大概分为两种:等价权函数法[5]和卡方检验法[6]。等价权函数法通过权函数来构造相应的等价权,对观测噪声方差阵进行降权处理,经常与自适应滤波[7-9]配合使用;卡方检验法的原理是服从高斯分布的状态量的马氏距离的平方服从卡方分布,通过设定阈值,对故障进行检测。在这两种抗差滤波方法上也有很多拓展:文献[10]中分析了由于观测值之间的相关性而导致的粗差误判问题,判断整体模型是否存在异常。文献[11]提出了错误预警率和故障探测率的概念,对故障的定位更加精确。文献[12]根据残差向量高斯分布统计特性对故障进行检测,对增益矩阵进行加权处理来抵消异常观测值影响。

上述抗差滤波算法都可以较好地解决观测值故障的问题,但也有可改进之处:两段判别域的故障检测方法中,被判定为无故障的区域解算出的导航信息也会出现偏差;当观测值出现超差情况时,当前时刻的观测值已经没有利用价值,简单的降权处理并不能保证观测值处在合理的范围内,应作废并构造一个新的观测值进行替换。据此,本文提出基于三段判断域的故障检测方法和最小二乘拟合的抗差滤波算法(以下简称拟合法),将检测区间分为三段:无故障、偏差和超差,当卫星导航观测值出现偏差时,对故障观测值进行降权处理;当检测到当前时刻卫星导航观测值有超差情况时,利用前n个时刻的卫星导航观测值构造的拟合函数进行一个时刻的外延,替换当前时刻的故障观测值。相对于传统的两段判别域的抗差滤波算法,三段判别域和基于最小二乘拟合的抗差滤波算法在导航精度和稳定性方面均有所提高。

1 以组合导航为模型的卡尔曼滤波算法

1.1 卡尔曼滤波

卡尔曼滤波[13]是20世纪60年代由数学家Kalman等人提出的基于线性最小方差估计的一种递推算法。

在许多的工程实践中,准确得到所需状态变量的真值是不可能的,只能通过观测信号和运动模型对状态变量进行估计和预测,来无限地接近真值,这也是卡尔曼滤波的基本思想[14-16]。卡尔曼滤波迭代流程如图1所示。

图1 经典卡尔曼滤波流程图

1.2 松组合模型

组合导航一般为惯性导航系统(inertial navigation system,INS)和全球卫星导航系统(global navigation satellite system,GNSS)[17-18]进行信息融合组成的导航系统[19-20],其组合方式大概分为3种:松组合、紧组合和深组合,基于算法层面的是松组合和紧组合。紧组合与松组合相比,状态量和观测量多了伪距和伪距率信息。松组合状态向量为

(1)

组合导航系统是一个复杂的系统,最简单的松组合系统也需要多信息融合才能运行,信号传递和数据处理过程中难免产生一些无用的噪声,所以导航系统一般使用卡尔曼滤波进行降噪处理[22]。以误差信息为状态量的松组合系统是一个线性系统,利用基本卡尔曼滤波便可以得到较好的效果。

2 基于卡方检验的故障检测方法

2.1 基于马氏距离的卡方检验法

卡尔曼滤波是一个迭代过程,残差向量包含有当前时刻全部的观测信息。残差向量表达式为

(2)

(3)

式中,μ取零向量;DM的平方服从卡方分布:

(4)

根据假设检验[24-26]的方法,设原假设为H0,假设导航系统观测值没有发生故障,也就是假设残差向量服从均值为零的高斯分布:

1.4 统计学方法 采用SPSS 22.0统计软件分析数据。计量资料以x±s表示，采用t检验；计数资料以百分比表示，采用χ2检验。以P<0.05为差异有统计学意义。

H0:DMk～χ(m)

(5)

式中，DMk为R时刻残差向量马氏距离的平方值。设显著性水平为α,即认为导航系统没有发生故障的概率为α,设

(6)

式中,χa(·)是显著水平为α的卡方分布。T(m)为检测量容许的最大值,也称阈值。当检测量DM

若当前时刻检测到故障时,可以将当前时刻的所有观测量当作故障值,但有可能会将其他正常的观测值当作故障值处理,造成的误差会在迭代的过程中逐渐积累,造成滤波发散。为防止这种现象的发生,本文采用单维度的故障检测方法[27],将设定的阈值变为

(7)

检测量的值变为

(8)

式中,i表示状态向量的第i维。

当检测量DM(i)

2.2 阈值的选取

DM(i)是根据当前时刻残差向量、估计误差方差阵等计算所得,而阈值T(1)需要根据卡方统计表选取,如表1所示。阈值选取的准确性决定了抗差滤波的准确性,表征对故障的容忍程度。阈值选取过大,含有较大粗差的观测值不容易被检测到,使用抗差滤波前后几乎无差别;阈值选取过小,一些相对正常的观测值会被当做故障观测值来处理,产生的误差经过迭代也会造成滤波发散。

如表1所示是显著性水平α和χ2(1)的分布表，从表中可以看出，显著性水平α越小，认为GNSS观测值发生故障的概率就越大，χ2(1)的值就越大，判断条件就越宽松。本文通过大量的实验，由检测量DM(i)的分布图分析得，检测量DM(i)的大小与观测噪声方差阵Rk和当前时刻粗差大小有关。如图2～图6所示为经度检测量分布图，这里仿真时长共500 s，设定每隔10 s加入粗差，提取这些时刻的卡方检测量的值作图。从图2～图6中可以看出，当误差小于5 m时，检测量数值有90%的概率小于2.706，当误差大于25 m时，检测量数值有96%的概率大于6.635。当误差为5～25 m时，检测量的值位于2.706～6.635之间，这个区间内有一段模糊区，检测量DM(i)穿插在卡方分布分位数的几个数值当中，仅用一个阈值判定故障与否很容易出现误判、漏判的情况。

图2 误差为5 m时卡方检测量分布图

图3 误差为10 m时卡方检测量分布图

图4 误差为15 m时卡方检测量分布图

图5 误差为20 m时卡方检测量分布图

图6 误差为25 m时卡方检测量分布图

3 三段判别域和最小二乘拟合的抗差滤波算法

拟合就是从一些数据点集(xi,yi)中找出规律，构造一个函数关系y=P(x)将这种规律表现出来，这个函数关系曲线要求尽可能靠近数据点，相比于插值[28]，拟合曲线不严格要求通过数据点[29-30]。Matlab软件提供了幂函数的最小二乘曲线拟合的功能，利用polyfit函数可以很快地求解出各个项的系数。

图7 三段判别域最小二乘拟合抗差滤波算法流程图

4 仿真分析

本文模拟一条飞机飞行轨迹进行仿真分析，飞行轨迹如图8所示，将飞行轨迹数据送入导航数据生成器，得到带有误差的模拟真实飞行环境的导航数据，仿真步长为1 s，总时长为501 s。惯性器件误差较小时的线性系统，一般仅利用位置观测值便可达到较好的效果，为简化模型，本文的仿真分析不使用速度观测值。

图8 模拟飞机飞行轨迹图

4.1 最小二乘拟合法抗差效果分析

拟合法是利用前n个时刻的GNSS观测值构造拟合曲线，进行一个时刻的外延，得到当前时刻的GNSS观测值。一般来说，n值越大，拟合效果越好，但计算量也会随之增大，本文经过大量实验分析得出n=4时拟合效果较好，计算量也不会很大。

为营造一个卫星信号丢失的环境，随机产生100个范围在0～501之间的整数作为观测值故障时刻，对这些时刻的GNSS东向、北向加入100 m的观测噪声，高度加入20 m观测噪声。

如图9和图10所示，一般的滤波方法不能消除观测值粗差带来的影响，在观测值发生故障的时刻，位置误差和速度误差都有较大波动。如图11和图12所示，INS在短期内的导航效果不错，误差较小，但由于INS积分误差会逐渐积累，随着时间的推移误差会越来越大。相比之下，使用最小二乘拟合的抗差滤波算法可以很好地减小故障观测值带来的影响，导航精度保持在一个良好的范围内。但拟合法抗差滤波也存在误差值高于INS的情况，这是由于基于拟合法的抗差滤波虽然可以修正故障观测值，但比较依赖于构造的函数模型，在某个时刻可能会出现较大误差。卡尔曼滤波是一个迭代的过程，每一次迭代都需要上个时刻的估计值作为输入，经过观测值的修正，曲线会逐渐趋于平稳。可以证明，最小二乘拟合的抗差滤波算法可以较为准确、稳定地输出导航信息。

图9 拟合法抗差滤波与无抗差滤波位置误差

图10 拟合法抗差滤波与无抗差滤波速度误差

图11 拟合法抗差滤波与纯惯导位置误差

图12 拟合法抗差滤波与纯惯导速度误差

4.2 三段判别域的故障检测方法仿真分析

当检测量的值小于2.706时，观测值误差小于5 m，认为观测值无故障，输出的导航信息也在可容忍的范围内。当检测量的值在2.706～6.635之间时，在两段判别域的故障检测方法中判定为无故障，在三段判别域的故障检测方法中判定为出现偏差，使用对观测值降权的方法来减弱故障的影响，以加入的观测值误差为20 m为例，如图13和图14所示。从图13和图14中可以看出，当观测值误差为20 m时，如不使用抗差滤波算法解算出的导航信息还是存在一定误差的，相对而言使用降权法抗差滤波时误差较小，具体误差范围如表2所示。当观测值误差为20 m时，检测量处在2.706～6.635这一区间的既有故障观测值，也有无故障观测值，但这些无故障观测值经过降权法抗差滤波之后的导航解算值仍处在合理范围内，所以在偏差情况使用降权法进行抗差滤波处理可以得到一个较好的结果。

图13 偏差情况降权法和无抗差滤波位置误差比较

图14 偏差情况降权法和无抗差滤波速度误差比较

抗差滤波无抗差滤波东向位置误差/m-11.8^6.06-25.8^32.6北向位置误差/m-6.64^7.54-20.5^34.9天向位置误差/m-2.62^2.57-7.9^13.9东向速度误/(m/s)-0.3^0.3-1.16^1.11东向速度误差/(m/s)-1.3^1.2-1.8^5.2东向速度误差/(m/s)-0.2^0.3-0.8^1.11

当检测量的值大于6.635时，观测值误差大于25 m，认为出现超差的情况，用基于最小二乘拟合的抗差滤波算法进行修正。当出现超差的情况时，理论上也可使用降权法对故障观测值进行修正。如图15和图16所示为每隔10个时刻观测值加入100 m的误差，拟合法与降权法的位置误差和速度误差图。从图15、图16及表3中可以看出，当故障时刻为单点时刻时，降权法和拟合法都可以减弱故障观测值带来的影响，由于拟合法构造拟合函数时的计算误差较大，降权法的抗差效果略好于拟合法，适用于丢包的情况。

表3 观测量误差为100 m时导航误差范围

但降权法也存在以下问题:当发生故障的时刻为连续的一段时间时,降权法的抗差效果并不理想。如图17和图18所示为10～200 s这一连续时间内观测值加入100 m的误差时,拟合法与降权法的位置误差和速度误差图。从图17和图18中可以看出,在10～200 s这一区间内,使用降权法进行抗差滤波,误差从0开始慢慢增大,随后误差值开始振荡,200 s后观测值正常,误差趋于平稳。相比之下,基于最小二乘拟合的抗差滤波算法在整个仿真区间位置误差和速度误差都比较稳定,抗差效果较好。归结原因主要为:拟合法在观测值发生故障时,将故障的GNSS观测值进行替换,使之处在相对合理的范围内,是从根源上解决观测值故障的一种方法,误差相对固定。降权法通过增大观测噪声方差阵Rk的值来减小滤波增益Kk的值,而Kk的值决定着估计均方误差阵Pk的大小,Pk经过迭代又作用于下一时刻的滤波增益值Kk+1,在连续时间的迭代中,会使增益阵渐渐失去合适的加权作用,造成发散。在单时刻故障中,当前时刻计算出的异常的Pk、Kk值会随着接下来几个时刻的正常值趋于正常,故导航效果良好。

图15 单时刻故障拟合法与降权法位置误差

图16 单时刻故障拟合法与降权法速度误差

图17 连续时刻故障拟合法与降权法位置误差

图18 连续时刻故障拟合法与降权法速度误差

5 结 论

针对导航系统在工作过程中GNSS观测值出现故障的问题,本文提出一种三段判别域的故障检测方法和基于最小二乘拟合的抗差滤波算法。三段判别域将故障分为无故障、偏差和超差3种情况,若判定观测值故障为偏差,用对观测值降权的方法进行抗差滤波处理,若故障为超差,用最小二乘拟合的方法进行抗差滤波。实验结果表明,三段判别域中对偏差情况进行降权处理比两段判别域中不做处理导航精度高。在超差情况下,单时刻故障时,使用降权法比拟合法精度略高,适用于丢包的情况。连续时间故障时使用拟合法比降权法导航精度高且更具稳定性,适用于链路中断的情况。