基于M-估计的抗欺骗卡尔曼滤波算法

张晓明，刘 磊

(尚禹河北电子科技股份有限公司，河北 石家庄 050091)

0 引言

GNSS作为国家重大基础设施的关键组成部分，已广泛应用于人民的日常生活中，例如民用航空、电力系统、时间同步网络、金融和交通运输等领域[1]。但是由于导航卫星信号的功率较低，易受各种有意或无意干扰的影响，而欺骗干扰由于具有功率较低、信号隐蔽性好的特点相对压制干扰对卫星导航系统的影响更大[2]，尤其是缺乏安全保护的民用导航信号[3]。

国内外的研究资料表明，关注GNSS欺骗干扰与反欺骗技术的学者越来越多，已成为GNSS领域内的研究热点[4-6]。目前常见的接收端欺骗干扰抑制技术主要有以下几种：① 通过对接收信号处理阶段信号的畸变来检测欺骗干扰[7-10]，然后剔除受欺骗影响的导航卫星，在多颗卫星被欺骗的情况下可能会导致参与定位的卫星数量不够；② 利用多个天线对比接收信号的差异或通过对欺骗干扰来波方向的估计检测和识别欺骗干扰[11-13]，比对信号差异的方法对运动的欺骗干扰信号无效，波达方向估计需要对每颗卫星信号相关后分别进行测向，测向后滤除或直接剔除受欺骗卫星，或再在欺骗方向形成零陷，计算量较大，不利于在接收机中实现；③ 组合导航的欺骗干扰检测和抑制方法[14-17]需要有惯导并且惯导和卫导必须为紧耦合或超紧耦合才能实现欺骗干扰抑制，而通用的小型接收机中一般没有惯导系统。因此，找到一种简单、易于在普通接收机中实现的欺骗干扰抑制方法成为目前的研究热点。

在此基础上，本文提出了一种应用在定位解算阶段的基于M-估计的抗欺骗扩展卡尔曼滤波算法，是对现在欺骗干扰检测与抑制算法的有益补充，解决了接收端的欺骗干扰抑制问题，又不增加现有接收机的硬件复杂度。

1 欺骗干扰

卫星导航的欺骗干扰具有和导航卫星信号相同的信号参数(信息码可能不同)，由于和真实信号功率相差不大，在干扰卫星导航接收机使其产生错误定位信息的同时还具备较好的隐蔽性。欺骗干扰主要有2种：产生式欺骗干扰和转发式欺骗干扰[18]，产生式欺骗干扰是根据侦察到的天上的卫星信号的码结构产生的与天上卫星信号相关性最大的伪随机码，然后调制上与导航电文结构一致的虚假电文。对于开放的民码信号，该种方式很容易实现，但是对于保密的军码信号，该种方式较难实现；转发式欺骗干扰是增加信号的传播时延，在干扰机收到卫星导航信号后，经过一定的滤波放大，再发射出去，该种欺骗干扰方式不需要区别军民信号，但是要注意收发之间的隔离。

欺骗干扰具有与导航信号相同的扩频码和载波频率的一种特殊信号，利用与真实信号的时延差或修改的导航电文影响导航接收机的正常定位，则欺骗干扰信号可表示为：

(1)

当存在欺骗干扰信号时，接收信号可表示为：

(2)

式中，n(t)为系统噪声；xi(t)为导航信号；L为当前可见卫星个数。

2 基于M-估计的抗欺骗卡尔曼滤波算法

若接收机捕获跟踪上欺骗信号，则将影响接收机的最终定位结果，因此，需要采用一定的滤波算法消除欺骗干扰对导航定位的影响。

2.1 M-估计在卫星导航定位中应用

(3)

(4)

从式(3)和式(4)可得，GNSS接收机定位解算的模型Zk=h[Xk,k]+Vk是非线性的，Vk为测量噪声，服从均值为零的高斯分布，k为时间序列。

将h[·]进行Taylor级数非线性展开，取一阶近似，则线性化的测量方程可表示为：

ΔZk=HΔXk+sk+vk,

(5)

式中，sk为欺骗干扰引起的误差；ΔZk为测量伪距，ΔXk为真实位置与估计位置之间的向量差；H是向量函数h[·]的雅克比矩阵，也称为观测矩阵，其元素为观测量对状态量的偏导。

残差f(·)函数和的线性化测量方程可写为：

(6)

式中，k时刻的测量量个数为n；f(·)为标量函数。

f(·)的最小化残差的和为：

(7)

进一步变形可得到：

(8)

式中，ψ(·)是f(·)对于相对状态矢量的导数，代表测量误差对定位结果的影响程度；hi是观测矩阵H=(hij)的第i行，即hi=hi1,hi2,…,him，m为参数的个数，进一步化简可得：

(9)

令：

(10)

则式(7)可化简为：

(11)

写为矩阵形式：

HTD(e)e=0，

(12)

式中，D(e)=diag[D(e1),D(e2),…，D(en)] 。

令e=HΔXk-ΔZk,则式(10)可写成：

(HTD(e)H)ΔXk=HTD(e)ΔZk。

(13)

式(13)与加权最小二乘方程一致，其解为：

(14)

为了使欺骗干扰对定位结果的影响最小，调整加权因子D(e)的值，使D(e)在测量误差越大时取值越小，不存在干扰时，所有权值取值为1。

考虑卫星导航接收机定位解算的特点，基于最小二乘的M-估计加权矩阵为：

(15)

2.2 结合M-估计的抗欺骗扩展卡尔曼滤波算法

将M-估计应用到卡尔曼滤波中，根据卡尔曼滤波的特点，基于卡尔曼滤波的M-估计加权矩阵可改写为：

(16)

式中，参数a和b为根据接收机实际应用场景设定的常数，数值的大小由接收机定位测速的误差决定。

扩展卡尔曼滤波的三维定位状态方程可写为：

GWk-1,

(17)

写成矩阵的形式为：

Xk=ΦXk-1+GWk-1，

(18)

则预测方程可表示为：

(19)

(20)

(21)

(22)

式中，Vk为测量余量；R为测量噪声的正定协方差矩阵；αk为测量余量协方差。

滤波增益为：

(23)

估计过程为：

(24)

(25)

3 计算机仿真与性能分析

仿真试验目的是为了验证M-估计的抗欺骗扩展卡尔曼滤波算法的抗欺骗性能。将本文提出的算法与最小二乘、扩展卡尔曼滤波和结合M-估计的抗欺骗扩展卡尔曼滤波算法在存在欺骗干扰时的定位性能对比。

仿真以 BDS系统为例，使用伪距单点定位方式，定位精度为米级，选取4颗GDOP值较好的卫星，测距误差为±3 m，假定状态矢量为(6 378 258,0,0,100,100,100,0,0)，初始滤波状态矢量为(6 378 000,4 000,2 000,0,0,0,0,0)。Q，R矩阵取值如下：

(26)

(27)

首先仿真欺骗干扰对伪距测量的影响，每间隔1 s采样一次，共计采样2 000 s，假设第500～860 s受到欺骗式干扰的影响，该欺骗干扰引起的伪距测量误差如图1所示，为1 200 m,大约相当于伪距偏移3个码片。

图1 欺骗引起的伪距误差Fig.1 Pseudorange error caused by spoofing

图2给出了第500～860 s受到欺骗干扰影响时，最小二乘算法的定位精度。

(a) X坐标轴引入的定位误差

(b) Y坐标轴引入的定位误差

(c) Z坐标轴引入的定位误差

从图2中可以看出,在受到欺骗干扰影响的第500～860 s，最小二乘算法在定位解算中不能消除欺骗干扰影响，欺骗干扰在X轴上引起的定位误差在2 000 m左右，在Y轴上引起的定位误差在-1 500 m左右，在高程上引起的定位误差在800 m左右。

图3给出了第500～860 s受到欺骗干扰影响时，传统扩展卡尔曼滤波算法的定位精度。

(a) X坐标轴引入的定位误差

(b) Y坐标轴引入的定位误差

(c) Z坐标轴引入的定位误差

从图3中可以看出，在受到欺骗干扰影响的第500～860 s，传统扩展Kalman滤波算法在定位解算中同样不能消除欺骗干扰影响。

图4给出了第500～860 s受到欺骗干扰影响时，结合M-估计的扩展卡尔曼滤波算法的定位精度。

(a) X坐标轴引入的定位误差

(b) Y坐标轴引入的定位误差

(c) Z坐标轴引入的定位误差图4 结合M-估计的扩展卡尔曼滤波法定位误差Fig.4 Positioning error of extended Kalman filtering combined with M-estimation

从图2～图4的对比中可以看出，最小二乘和扩展卡尔曼滤波算法都不能抑制欺骗干扰对导航定位的影响，定位结果出现了较大的误差，将近2 000 m。本文提出算法在存在欺骗干扰的条件下，定位误差优于10 m，说明本文算法对欺骗干扰具有较好的抑制性能，定位精度较好。

4 结束语

本文提出了一种基于M-估计的抗欺骗扩展卡尔曼滤波算法，结合稳健统计中的M-估计来调整扩展卡尔曼滤波的状态更新过程，消除突发欺骗干扰对导航定位的影响，提升接收机定位解算过程中的欺骗干扰抑制性能。通过与最小二乘定位和统扩展卡尔曼滤波定位算法的对比表明，在引入1 200 m的伪距误差时，最小二乘定位和统扩展卡尔曼滤波定位算法在3个方向上均出现了较大的定位误差，而本文提出的抗欺骗干扰定位算法在3个方向定位精度均优于10 m。本算法适用于对成本敏感且具有一定抗欺骗干扰能力的普通接收机，可以在不增加硬件成本的前提下使接收机具有欺骗干扰抑制能力。