1.   引言

全球导航卫星系统(Global Navigation Satellite System, GNSS)的定位精度高和实时性好等优点使其在民航航行新技术中扮演着重要角色^[1]。然而，GNSS信号功率低和民用粗码结构公开等因素导致其极易受到欺骗式干扰的影响^[2,3]。欺骗式干扰使用假冒的卫星信号，致使机载接收机生成错误的定位结果，这严重影响了飞机的飞行安全。软件无线电技术(Software Defined Radio, SDR)与硬件平台的发展使得个人生成模拟GNSS信号将更加容易^[4,5]，未来蓄意破坏者对目标用户进行有意干扰的频率将大大增加。GNSS欺骗逐渐成为不容忽视的威胁。

针对欺骗式干扰，主要有信号功率监视^[6]、密码认证^[7-9]、天线阵列处理^[10,11]和使用辅助导航信息^[12]等检测方式。惯性导航系统(Inertial Navigation System, INS)是机载自主式导航设备，不受无线电信号的干扰，与GNSS有较高的互补性，是检测机载GNSS欺骗最常用的辅助导航工具^[13,14]。基于卡尔曼滤波新息序列的欺骗检测算法，不需要对基础设施进行修改，也不需要增加额外的设备，利用飞机上已有的INS/GNSS数据即可实现检测，便于飞行过程中使用^[15]。欺骗式干扰会导致INS/GNSS系统量测值发生幅值较大的阶跃变化，此时新息序列的均值不为0，通过该新息序列构造的卡方检验统计量会超过阈值，从而检测到欺骗式干扰的存在^[16]。也就是说，卡方检验算法有意义的条件是构成卡方检验统计量的新息序列的均值能够在存在欺骗式干扰的情况下始终不为0。然而，当量测值突变幅度保持不变时，卡尔曼滤波的估计值和一步预测值会由于欺骗量测值的修正而逐渐和量测值“一致”起来，新息序列均值重新归0，检验统计量重新回到阈值之下，卡方检测算法失效，即使后面存在欺骗，检测算法也无法通过统计量与阈值的关系进行判断。因此，INS/GNSS系统无法依据传统卡方检测算法来判断给出的定位信息是否正确。另外，以上使用辅助导航信息进行欺骗检测的方法大多局限于INS/GNSS组合方式，并没有考虑到机载多导航源这一特点。

机载测距机(Distance Measuring Equipment, DME)是一种由地面应答机和机载询问机两部分组成的无线电测距导航系统。该系统通过测量脉冲信号在两者之间的传播时间来获取斜距信息，进而确定飞机的3维位置^[17]。GNSS欺骗信号的功率比真实GNSS信号功率略高且在接收机热噪声下，DME信号的功率远高于接收机热噪声，因此DME系统不会受到欺骗式干扰的影响^[18]。本文使用联邦卡尔曼滤波进行INS,GNSS和DME的数据融合，引入融合后的INS/DME数据参与构建欺骗检验统计量。当机载GNSS欺骗达到一定程度时，本文算法能够给出较为准确的欺骗持续时间。本文进一步将所设计的算法应用于组合导航系统中，得到了正确的定位信息，最后通过仿真实验对所提算法的有效性进行了验证。

2.   组合导航模型与问题描述

本节分别给出了松组合方式下INS/GNSS系统和INS/DME系统的卡尔曼滤波的状态方程和量测方程，并对欺骗式干扰环境下的INS/GNSS量测方程进行量化分析，简要说明了传统检测算法的失效原因。

2.1   状态方程与量测方程的构建

使用导航系原点固定的INS直接输出北、东和下3个方向的位置信息，以INS的误差参数作为状态量，得到卡尔曼滤波的状态方程为

其中，${\phi _{k/k - 1}}$表示从$k - 1$到$k$时刻的状态转移矩阵；${{\boldsymbol{\varGamma}} _{k - 1}}$和${{\boldsymbol{w}}_{k - 1}}$分别表示$k - 1$时刻的噪声分配矩阵和过程噪声矩阵；状态量${{\boldsymbol{X}}_k}$是$15 \times 1$的矩阵，包括

其中，$\delta {\boldsymbol{P}}_k^{({\text{n}})}$和$\delta {\boldsymbol{v}}_k^{({\text{n}})}$分别表示INS的位置累积误差和速度累积误差；${\boldsymbol{\varphi}}_k^{({\text{n}})}$表示INS的失准角误差；$\delta {\boldsymbol{f}}_k^{({\text{b}})}$表示载体坐标系下的比力误差矢量；$\delta {\boldsymbol{w}}_{{\text{ib,}}k}^{({\text{b}})}$为载体坐标系下载体坐标系相对于惯性坐标系的旋转角速度的误差矢量。

假设飞机的真实位置为${{\boldsymbol{P}}_k}$, INS, GNSS和DME的定位结果分别为

其中，$\delta {{\boldsymbol{P}}_k}$表示INS位置的累积误差；${\boldsymbol{n}}_k^{{\text{GNSS}}}$和${\boldsymbol{n}}_k^{{\text{DME}}}$分别表示GNSS和DME的位置解算误差，物理含义为0均值的高斯白噪声。

式(3)减去式(4)，得到卡尔曼滤波量测方程1为

其中，${{\boldsymbol{H}}_k}$为量测系数矩阵，用于提取状态量中的位置误差$\delta {{\boldsymbol{P}}_k}$，表示为

其中，${{\boldsymbol{I}}_{3 \times 3}}$为$3 \times 3$的单位矩阵；${{{ {\textit{0}}}}_{3 \times 12}}$为$3 \times 12$的0矩阵。式(3)减去式(5)，得到卡尔曼滤波量测方程2为

根据式(1)与式(6)，得到INS/GNSS 组合下卡尔曼滤波的一步预测值和估计值分别为${\boldsymbol{\hat X}}_{k/k - 1}^{{\text{INS}}/{\text{GNSS}}}$和${\boldsymbol{\hat X}}_k^{{\text{INS/GNSS}}}$。同理，根据式(1)与式(8)，得到INS/DME组合下卡尔曼滤波的一步预测值和估计值分别为${\boldsymbol{\hat X}}_{k/k - 1}^{{\text{INS/DME}}}$和${\boldsymbol{\hat X}}_k^{{\text{INS/DME}}}$。

2.2   INS/GNSS量测方程分析

假设目标接收机在$k \in \left[ {t,t + n} \right]$受到欺骗，这段时间内的INS/GNSS量测值比正常情况下多一个偏差$\delta {\boldsymbol{P}}_k^{{\text{bias}}}$，表示为

刚加入欺骗时，由于$\delta {\boldsymbol{P}}_k^{{\text{bias}}}$的存在导致新息序列不满足0均值这一统计特性，使用新息序列构造的瞬时卡方检验统计量超过阈值，检测到欺骗存在。然而，${\boldsymbol{\hat X}}_{k/k - 1}^{{\text{INS}}/{\text{GNSS}}}$和${\boldsymbol{\hat X}}_k^{{\text{INS/GNSS}}}$在卡尔曼滤波的循环迭代中不断受到${\boldsymbol{Z}}_k^1$修正，导致其不再是对${{\boldsymbol{X}}_k}$的一步预测和估计，而是对$\left( {{{\boldsymbol{X}}_k} + {\boldsymbol{X}}_k^{{\text{bias}}}} \right)$ 整体进行一步预测和估计。新息序列重新满足0均值这一统计规律，即使后面仍然存在欺骗，瞬时检验统计量也会回到阈值之下。因此，传统算法会出现无法根据统计量与阈值的关系来判断有无欺骗的情况。

根据上述分析可知，传统算法是否有效关键在于新息序列的统计特性是否发生变化。而新息序列的统计特性是否变化则取决于量测值是否发生阶跃变化。当新息序列重新满足0均值后，$\delta {\boldsymbol{P}}_k^{{\text{bias}}}$突然减小(其中$\delta {\boldsymbol{P}}_k^{{\text{bias}}}$突然减小到0意味着欺骗结束，是一种特殊情况)或增大都会导致量测值再次发生突变，而此时的统计量会再次超过阈值。因此，当统计量第2次超过阈值后，不一定就是由于欺骗结束而导致，进而说明，传统算法无法依据统计量超过阈值的次数来判断欺骗是否结束。也就是说，传统算法存在检测不到欺骗持续时间的局限性。

3.   欺骗检测算法

本文主要由传感器模块、欺骗检测模块和定位解算模块构成，如图1所示。传感器模块提供检测所必需的INS, GNSS和DME位置信息。欺骗检测模块给出了基于重构新息序列的有限记忆卡方检测算法的构造过程并证明了该方法的可行性。在此基础上，本文进一步利用DME数据，给出定位解算模块。该模块能够在存在GNSS欺骗的情况下提供正确的飞行引导。下面对新的检测算法进行详细介绍。

图 1  新的检测算法原理图

下载: 全尺寸图片 幻灯片

3.1   重构新息序列

INS/DME系统的${\boldsymbol{Z}}_k^2$不受GNSS欺骗的影响，因此${\boldsymbol{\hat X}}_{k/k - 1}^{{\text{INS}}/{\text{DME}}}$和${\boldsymbol{\hat X}}_k^{{\text{INS}}/{\text{DME}}}$始终是对${{\boldsymbol{X}}_k}$进行一步预测和估计。利用已有${\boldsymbol{Z}}_k^1$和${\boldsymbol{Z}}_k^2$重新构造一个量测值${\boldsymbol{Z}}_k^{{\text{new}}}$，并与${{\boldsymbol{H}}_k}{\boldsymbol{\hat X}}_{k/k - 1}^{{\text{INS}}/{\text{DME}}}$做差构建新息序列，该新息序列不参与卡尔曼滤波过程，满足了在有欺骗下新息序列均值不为0的要求。如果构造的${\boldsymbol{Z}}_k^{{\text{new}}}$与${\boldsymbol{Z}}_k^2$的协方差矩阵相同，则与${{\boldsymbol{H}}_k}{\boldsymbol{\hat X}}_{k/k - 1}^{{\text{INS}}/{\text{DME}}}$做差时更加一致，并且在后面构建检验统计量时可以直接使用已有的INS/DME的新息协方差矩阵${\boldsymbol{C}}_k^{{\text{DME}}}$作为${\boldsymbol{C}}_k^{{\text{new}}}$。

基于以上分析，使用$3 \times 3$的系数矩阵${{\boldsymbol{A}}_1}$和${{\boldsymbol{A}}_2}$，构造新的量测值${\boldsymbol{Z}}_k^{{\text{new}}}$为

在有、无欺骗的不同时间段内，${\boldsymbol{Z}}_k^{{\text{new}}}$可以用分段函数表示为

令${\boldsymbol{n}}_k^{{\text{new}}} = {{\boldsymbol{A}}_1}{\boldsymbol{n}}_k^{{\text{GNSS}}} + {{\boldsymbol{A}}_2}{\boldsymbol{n}}_k^{{\text{DME}}}$，并进行如式(12)的约束

其中，${\boldsymbol{R}}_k^{{\text{GNSS}}}$, ${\boldsymbol{R}}_k^{{\text{DME}}}$和${\boldsymbol{R}}_k^{{\text{new}}}$分别为${\boldsymbol{n}}_k^{{\text{GNSS}}}$, ${\boldsymbol{n}}_k^{{\text{DME}}}$和${\boldsymbol{n}}_k^{{\text{new}}}$的协方差矩阵。通过式(12)解算出${{\boldsymbol{A}}_1}$和${{\boldsymbol{A}}_2}$并代回式(10)即可得到已知的新的量测值。根据该量测值，重新构造新息序列$\delta {\boldsymbol{z}}_k^{{\text{new}}}$为

在式(13)的基础上进一步构造检验统计量，用于实现对欺骗开始和结束时间的较为准确的获取。下面给出有限记忆检验统计量的具体构造过程和相关的判决准则。

3.2   构建有限记忆卡方检验统计量

使用新息序列可以构造成瞬时检验统计量和累加检验统计量^[19]。根据式(13)构建瞬时检验统计量虽然能够对GNSS欺骗起到及时报警的作用，但由于每一时刻的统计量取值不尽相同，当${{\boldsymbol{A}}_1}\delta {\boldsymbol{P}}_k^{{\text{bias}}}$较小时会出现检验统计量在阈值上下波动的情况，无法获得较为准确的欺骗持续时间。根据式(13)构造累加检验统计量可以避免波动的问题。然而，累加检验统计量包含了大量老旧的数据，无法及时反映当前样本的分布规律，导致欺骗结束后的统计量无法及时回到阈值之下，因此使用累加检验统计量得到的欺骗持续时间远远大于真实的欺骗持续时间。

为此，本文利用滑动窗口算法构造有限记忆卡方检验统计量。有限个瞬时检验统计量的累加可以避免统计量在阈值上下波动的情况，同时能够保证在欺骗结束之后统计量快速回到阈值之下。

有限记忆卡方检验统计量$s_k^2$的表达式为

其中，$N$表示滑动窗口的长度；协方差阵${\mathbf{C}}_i^{{\text{new}}} = {{\mathbf{H}}_i}{\mathbf{P}}_{i/i - 1}^{{\text{DME}}}{\mathbf{H}}_i^{\text{T}} + {{\mathbf{R}}^{{\text{new}}}} = {\mathbf{C}}_i^{{\text{DME}}}$。在没有欺骗的情况下，$s_k^2$符合如式(15)的分布规律

其中，$n$为自由度，表示为

在给定恒虚警率${P_{{\text{FA}}}}$的情况下，阈值可以通过卡方分布的反函数确定，即

其中，卡方分布函数$F(x|n)$表达式为

其中，$\varGamma \left( x \right)$是gamma函数。根据已知的阈值$T_k^2$和检验统计量$s_k^2$可以进行如下判决。

如果存在GNSS欺骗，则有

反之，则不存在欺骗式干扰。使用基于重构新息序列的有限记忆卡方检验统计量能够较为准确地判断欺骗的开始和结束，为使用GNSS还是DME进行准确的飞行引导提供了判断依据。

3.3   新的检测算法有效性证明

根据2.2节的分析可知，新息序列的均值是否为0，直接决定了后续构造的检验统计量能否检测到欺骗。下面给出新的检测算法的有效性证明。

将式(12)代入式(11)，对新的量测值进行化简，得到

其中，${{\boldsymbol{A}}_1}$和${{\boldsymbol{A}}_2}$是通过式(12)计算后的矩阵。将式(20)代入式(13)，得到不同时间段的重构新息序列为

令${\boldsymbol{\tilde X}}_{k/k - 1}^{{\text{INS/DME}}} = {{\boldsymbol{X}}_k} - {\boldsymbol{\hat X}}_{k/k - 1}^{{\text{INS/DME}}}$，表示状态估计的误差。将式(21)进一步化简得到

一步预测值${\boldsymbol{\hat X}}_{k/k - 1}^{{\text{INS/DME}}}$始终是对${{\boldsymbol{X}}_k}$的无偏估计，即$E\left( {{\boldsymbol{\hat X}}_{k/k - 1}^{{\text{INS/DME}}}} \right){\text{ = }}{{\boldsymbol{X}}_k}$。状态估计误差${\boldsymbol{\tilde X}}_{k/k - 1}^{{\text{INS/DME}}}$的均值可化简为

进一步得到重构新息序列的均值表达式为

当存在欺骗时，$\delta {\boldsymbol{P}}_k^{{\text{bias}}} \ne {{{ {\textit{0}}}}_{3 \times 1}}$，则始终有重构新息序列的均值$E\left( {\delta {\boldsymbol{z}}_k^{{\text{new}}}} \right) = {{\boldsymbol{A}}_1}\delta {\boldsymbol{P}}_k^{{\text{bias}}} \ne {{{ {\textit{0}}}}_{3 \times 1}}$，保证了在有欺骗情况下算法的有效性；当欺骗结束时，有$\delta {\boldsymbol{P}}_k^{{\text{bias}}} = {{{ {\textit{0}}}}_{3 \times 1}}$，新息序列均值重新为0，保证了在欺骗结束后检验统计量能够回到阈值之下。

4.   定位算法

当第3节所提算法的检验统计量不高于阈值时，将INS/GNSS和INS/DME进行数据融合，从而得到最优估计状态的协方差矩阵和最优估计值^[20]。

最优估计值的协方差矩阵${\boldsymbol{P}}_{k/k}^{{\text{INS}}/{\text{GNSS/DME}}}$表示为

其中，${\boldsymbol{P}}_{k{\text{/}}k}^{{\text{INS/GNSS}}}$和${\boldsymbol{P}}_{k{\text{/}}k}^{{\text{INS/DME}}}$分别为INS/GNSS和INS/DME系统的估计值均方误差矩阵。

在式(25)的基础上，得到系统的最优估计值表达式为

当检验统计量高于阈值时，对GNSS进行系统隔离，只使用INS/DME对INS的累积误差进行估计。由此得到有、无GNSS欺骗下，组合导航系统的位置解${\boldsymbol{P}}_k^{{{{\rm{correct}} - {\rm{INS}}}}}$为

5.   仿真实验

本节通过仿真来验证新算法的有效性。模拟一段时间为180 s的直线飞行轨迹，初始位置为北纬40°，东经116°，高度1000 m。该实验使用GNSS模拟器生成GNSS信号，然后用软件接收机接收得到GNSS的位置信息。按照实际应用中对INS传感器精度的要求，设置陀螺仪随机噪声和加速度计随机噪声分别为0.07 °/h²和3 mg，设置陀螺仪和加速度计的1阶马尔可夫过程噪声分别为0.08 °/h²和12 mg，设置马尔可夫过程的相关时间为3600 s。根据真实DME测量精度，设置其测量误差为180 m。

本文设置欺骗的开始时间和结束时间分别为第40 s和第141 s。默认使用的欺骗信号在北、东、地3个方向上引起的偏差相同，以东方向位置误差为例，如图2所示。由图2可知，如果接收机接收真实GNSS信号，则接收机解算出来的位置误差均值始终为0。如果软件接收机从第40 s开始接收欺骗信号，可以看到，在第40～140 s这一时间段内，接收机解算的位置误差均值为300 m。由此说明，接收机接收欺骗信号解算的位置与接收真实GNSS信号解算位置存在300 m的偏差。

图 2  接收机有、无欺骗的东向位置误差

下载: 全尺寸图片 幻灯片

使用传统瞬时卡方算法对图2的欺骗场景进行检测，如图3所示。由图3可以看到，传统算法在第40 s超过阈值而报警，但在第53 s之后，传统算法的检验统计量重新回到阈值之下。第53～140 s这一时间段内，虽然存在GNSS欺骗但不会出现警报，直到第141 s欺骗结束引起量测值第2次发生突变，检验统计量再次超过阈值而告警。由此不难得出，传统算法给出的统计量与阈值的关系不能用来判断欺骗是否结束。

图 3  基于传统新息序列的瞬时卡方检测算法

下载: 全尺寸图片 幻灯片

在相同的实验条件下，使用本文设计的算法对图2的欺骗场景进行检测，如图4所示。滑动窗口的长度$N$设置为10，有限记忆卡方检验统计量在第41 s超过阈值而发出警报，表明该算法检测到了欺骗式干扰的存在。在第149 s检验统计量重新回到阈值之下，表明本文算法检测到了GNSS欺骗结束。由此可知，新的检测算法可以依靠统计量与阈值的关系来判断欺骗的开始和结束。

图 4  基于重构新息序列的有限记忆卡方检测算法

下载: 全尺寸图片 幻灯片

图5分别给出了INS/GNSS系统位置误差、INS/DME系统定位误差和本文算法的误差曲线，其中3个子图分别表示北、东、地3个方向。由图可知，在没有欺骗式干扰时，INS/GNSS系统的误差曲线均值为0。第40～141 s这一时间段内，GNSS欺骗导致INS/GNSS系统的误差曲线被带到0均值以外。无论是否存在GNSS欺骗，INS/DME系统误差曲线始终保持0均值。然而，DME系统本身的定位误差远远大于GNSS，从而导致INS/DME的定位误差大于INS/GNSS。当本文所提的检测算法没有超过阈值时，使用INS/GNSS/DME作为组合导航系统的定位解。当检测到GNSS欺骗后，本文的定位算法立刻停止GNSS参与系统定位，因此红色的误差曲线能够始终保持0均值的统计规律。在GNSS系统隔离期间，本文的定位算法实质上就是INS/DME，组合导航系统的定位精度会发生下降。然而，无论是否存在欺骗式干扰，本文给出的定位方案都优于单一使用INS/GNSS或INS/DME。

图 5  组合导航系统位置误差

下载: 全尺寸图片 幻灯片

下面分析不同DME测量精度下，欺骗引起的位置偏差对新算法检测概率的影响。接收机分别接收不同的欺骗信号，使其解算位置分别偏离真实GNSS位置20 m, 25 m, 30 m直到100 m。在存在GNSS欺骗的时间段内，如果出现检验统计量超过阈值的情况，则视该次实验检测到GNSS欺骗。进行1000次蒙特卡罗仿真，统计传统算法和本文算法检测到欺骗的次数，如图6所示。由图6可知，DME的定位精度越高，能够检测到的位置偏差量越小，检测概率越高。如果DME的测量误差为20 m，当3个方向都偏离20 m时，本文算法的检测概率就达到100%。相同实验条件下，传统算法需要当3个方向都偏离70 m时检测概率才能达到100%。由此说明，如果未来DME的定位精度提高或者存在其他定位精度更高的辅助导航源，所提算法都将有比现在更好的检测效果。

图 6  不同DME测量误差下本文算法的检测概率

下载: 全尺寸图片 幻灯片

上面实验提到，新算法的滑动窗口长度设置为10，下面给出$N{\text{ = 10}}$的选取依据。$N$的选取需要综合考虑欺骗引起的位置偏差量、漏警概率和时间延迟量3个因素。接收机分别接收不同的欺骗信号，使其解算位置分别偏离真实GNSS位置200 m, 250 m, 300 m直到500 m。在存在GNSS欺骗的时间段内，如果检验统计量出现低于阈值的情况，则视该次实验出现漏警。进行1000次蒙特卡罗实验，对位置偏差量和滑动窗口长度共同影响下的新算法的漏警概率进行统计，如图7所示。由图7可知，当本文算法的滑动窗口长度$N$设置为1时，其漏警概率在位置偏差为500 m时才能达到0%。相同实验条件下，当本文算法的滑动窗口长度$N$设置为10时，其漏警概率在加入250 m的GNSS欺骗时就能够达到0%。对图7进一步分析可知，如果GNSS欺骗为250 m及以上时，新算法的$N \ge10$就可以保证漏警概率为0%。而对于200 m及以下的GNSS欺骗，需满足新算法的$N \ge {\text{20}}$才可以达到0%的漏警概率。

图 7  不同滑窗长度下本文算法的漏警概率

下载: 全尺寸图片 幻灯片

进一步分析不同程度的欺骗下，滑动窗口长度对新算法检测延迟的影响，如图8所示。由图8可知，当$N \ne {\text{1}}$时，滑动窗口的存在会造成本文算法在欺骗检测上出现延迟。对于任意一个GNSS欺骗，滑动窗口的长度$N$越大，两幅子图的时间延迟量越大，所设计的检测算法给出的欺骗持续时间越不准确。假如欺骗式干扰至少会引起250 m的位置偏差，在满足漏警概率为0%的条件下，应让$N$尽量小。综合欺骗引起的位置偏差量、漏警概率和时间延迟量3个因素，本文选取滑动窗口的长度$N$为10。

图 8  不同滑窗长度下本文算法的检测延迟

下载: 全尺寸图片 幻灯片

6.   结束语

基于传统新息序列的卡方检测算法无法给出欺骗的持续时间，从而导致INS/GNSS系统无法根据该算法的检测结果判断其定位解算是否正确。本文提出一种基于重构新息序列的有限记忆卡方检测算法。本算法使用INS,GNSS和DME的数据构造一种不参与卡尔曼滤波的新息序列，然后将该新息序列构造成有限记忆卡方检验统计量，从而实现对欺骗式干扰的检测。仿真表明，当欺骗式干扰造成250 m及以上的位置偏差时，所设计的算法能够实现对欺骗持续时间较为准确的获取。本文进一步将所提算法应用于INS/GNSS/DME系统中，给出了正确的飞行引导方案。由检测概率的仿真图可知，如果未来DME的定位精度提高，或者存在其他定位精度更高的导航源，本文算法都将有比现在更好的检测效果和定位效果。