一种组合导航系统的欺骗干扰检测方法与装置

摘要

本发明提出一种组合导航系统的欺骗干扰检测方法与装置，属于卫星导航技术领域。所述方法包括：构造零速检测统计量；构造滤波模型；构造静止状态欺骗检测统计量；构造运动状态欺骗检测统计量。本发明将抖动的欺骗速率建模为服从非零均值高斯分布的随机变量，通过比较静止无欺骗状态下统计的参考噪声方差和运动状态下的新息序列方差，构建新息序列方差的检测统计量，进行欺骗式干扰检测，对欺骗速率均值不敏感，针对缓变式欺骗干扰具有很好的检测性能。(56)对比文件史密;陈树新;吴昊;毛虎.拒止环境实现注入的GPS欺骗干扰.空军工程大学学报(自然科学版).2015,(06),27-31.

权利要求

1.一种组合导航系统的欺骗干扰检测方法，其特征在于，所述组合导航系统包括全球卫星导航系统GNSS接收机和惯性导航系统INS，所述惯性导航系统INS包括测量惯性单元IMU；所述方法包括：

步骤S1、获取所述测量惯性单元IMU输出的三维加速度和角速度测量值，作为检测统计量，并基于所述检测统计量构造零速检测统计量，所述零速检测统计量用于确定载体是否处于静止状态；

步骤S2、构造滤波模型，所述滤波模型用于估计当前系统状态，并构造欺骗检测统计量，其中，所述滤波模型为卡尔曼滤波模型；

步骤S3、响应于所述载体处于所述静止状态，根据所述当前系统状态构造所述静止状态下的欺骗检测统计量，所述静止状态下的欺骗检测统计量用于确定所述载体在所述静止状态下是否受到欺骗干扰，并且响应于所述载体在所述静止状态下未受到所述欺骗干扰，进一步确定参考噪声；

步骤S4、计算所述参考噪声的方差，利用所述参考噪声的方差构造运动状态下的欺骗干扰统计量，所述运动状态下的欺骗干扰统计量用于判断所述载体在所述运动状态下是否受到欺骗干扰；

其中，所述载体为车辆。

2.根据权利要求1所述的一种组合导航系统的欺骗干扰检测方法，其特征在于：在所述步骤S1中：

利用所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值，作为所述检测统计量，构造广义似然比公式，得到零速检测公式：其中， 分别为所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值的测量方差，Zn为第n次检测的检测统计量集合，所述检测统计量为所述测量惯性单元IMU在检测时间段内各时刻的测量值，具体包括第k测量时刻的三轴加速度 和角速度N为检测窗口内数据个数，g为重力常量， 为检测窗口长度内的加速度平均值， 为加速度平均值 的幅值，T表示矩阵转置，Ωn为检测窗口长度， 表示

的二范数运算， 表示 的二范数运算；

其中，当所述检测统计量T(zn)小于设定的门限γ时，确定所述载体处于所述静止状态；

在所述步骤S2中：

所述组合导航系统采用所述卡尔曼滤波模型，第k个测量时刻的状态向量为：T

xk＝[δpk，δvk，δφk，δωk，δfk，bclk] 其中， 表示位置误差向量，由经度、纬度和高度组成；δvk＝[δvE，k，δvN，k，δvC，k]表示地球参考速度误差向量；δφk＝[δpk，δrk，δAk]为姿态误差向量；δωk＝[δωR，k，δωF，k，δωU，k]为陀螺仪零偏误差向量；δfk＝[δfR，k，δfF，k，δfU，k]为加速度计零偏误差向量；bclk＝[δbr，δdr]为接收机钟差和钟漂；其中，下标E，N，C分别表示当地地理坐标系中的东、北、天分量；下标R，F，U分别表示载体坐标系中的右、前、上分量；第k个测量时刻的测量向量zk表示如下：其中， 分别为k时刻第i颗卫星的观测伪距和观测伪距率，i＝1，2，L，M，M为观测卫星数量；则系统状态方程和测量方程为：其中，xk和xk‑1分别表示k和k‑1时刻的状态向量，wk‑1表示k‑1时刻的状态噪声向量，Φk/k‑1为k‑1时刻的状态转移矩阵，由所述惯性导航系统INS的误差方程和所述全球卫星导航系统GNSS接收机的钟差模型得出，连续到离散过程采用一阶近似；Hk表示k时刻测量矩阵，vk为k时刻观测噪声向量，wk为k时刻状态噪声向量，wk和vk皆假设为零均值的高斯白噪声，协方差分别用Qk和Rk表示；

在卡尔曼滤波的基础上，利用当前时刻的测量向量zk和状态预测向量 构建新息向量δzk为：

3.根据权利要求2所述的一种组合导航系统的欺骗干扰检测方法，其特征在于，在所述步骤S3中：

当检测到所述载体处于所述静止状态时，所述载体的车体坐标系b系下三个方向的速b

度为零，v＝0，在该k时刻由所述全球卫星导航系统GNSS接收机解算得到的所述载体的速度测量值为 对所述全球卫星导航系统GNSS接收机进行欺骗干扰检测，利用所述全球卫星导航系统GNSS接收机解算得到的所述载体的速度测量值作为所述静止状态下的欺骗检测统计量，构造广义似然比公式，得到所述静止状态下的欺骗干扰检测公式：其中，zn为第n次检测的检测统计量集合，N为检测窗口内数据个数，Ωn为检测窗口长度， 为GNSS观测方差， 表示的二范数运算；

当所述静止状态下的欺骗检测统计量T(zn)小于设定的门限γ时，确定所述全球卫星导航系统GNSS接收机未受到欺骗干扰，获取未受欺骗干扰时段的新息序列δzk作为所述参考噪声，用于所述运动状态下的欺骗干扰检测，所述未受欺骗干扰时段的新息序列的均值为零，统计不同卫星新息序列的方差为：其中， 为第j颗卫星新息的方差，M为统计新息个数，ΩM为未受欺骗干扰的窗口长度，δzk，j为k时刻第j颗卫星的新息。

4.根据权利要求3所述的一种组合导航系统的欺骗干扰检测方法，其特征在于，在所述步骤S4中：

在所述静止状态下得确定所述参考噪声方差，并且当所述载体处于所述运动状态时，将欺骗速率作为满足非零均值高斯分布的随机变量，设所述运动状态下k时刻第j颗可见卫星的伪距率新息为 当所述全球卫星导航系统GNSS接收机未受到所述欺骗干扰时，伪距率新息方差应为所述静止状态下统计的对应参考噪声方差，当所述全球卫星导航系统GNSS接收机受到所述欺骗干扰时，欺骗速率的抖动使得新息方差大于参考噪声方差；令H0为无欺骗状态，H1为欺骗状态，则所述运动状态下的欺骗干扰检测为二元假设模型：

其中，检测序列长度为M，w[k]为已知参考噪声方差 的零均值高斯噪声，s[k]为由欺骗信号引入的方差 未知的零均值高斯分布信号；则广义似然比为：采用MLE方法估计 通过取对数得到：对 求导得到：

令导数为0，得到欺骗信号方差的MLE估计为：若 则相应的MLE为 与参数约束保持一致，则 的MLE为令 其中上标+表示

为正时，其为MLE；

当 时，判定为H1，将估计值 代入对数似然比可得：

‑1

令 g(z)＝z‑lnz‑1为单调递增函数，对于z＞1，逆g 存在，若或 判定为H1，欺骗检测统计量

为对欺骗信号方差的MLE估计值

2

在H0条件下，新息序列满足零均值方差为σj的高斯分布，当样本数M足够时，新息方差渐进满足高斯分布：

其中，

则 的概率密度函数为：

根据 的概率密度函数和虚警概率计算门限，门限γ为：‑1

其中，Q 为逆高斯右尾概率函数，PFA为虚警概率，当统计量 大于设定的门限γ时判定所述载体的所述全球卫星导航系统GNSS接收机受到欺骗干扰。

5.一种组合导航系统的欺骗干扰检测装置，其特征在于，所述组合导航系统包括全球卫星导航系统GNSS接收机和惯性导航系统INS，所述惯性导航系统INS包括测量惯性单元IMU；所述装置包括：

第一处理单元，被配置为：获取所述测量惯性单元IMU输出的三维加速度和角速度测量值，作为检测统计量，并基于所述检测统计量构造零速检测统计量，所述零速检测统计量用于确定载体是否处于静止状态；

第二处理单元，被配置为：构造滤波模型，所述滤波模型用于估计当前系统状态，并构造欺骗检测统计量，其中，所述滤波模型为卡尔曼滤波模型；

第三处理单元，被配置为：响应于所述载体处于所述静止状态，根据所述当前系统状态构造所述静止状态下的欺骗检测统计量，所述静止状态下的欺骗检测统计量用于确定所述载体在所述静止状态下是否受到欺骗干扰，并且响应于所述载体在所述静止状态下未受到所述欺骗干扰，进一步确定参考噪声；

第四处理单元，被配置为：计算所述参考噪声的方差，利用所述参考噪声的方差构造运动状态下的欺骗干扰统计量，所述运动状态下的欺骗干扰统计量用于判断所述载体在所述运动状态下是否受到欺骗干扰；

其中，所述载体为车辆。

6.根据权利要求5所述的一种组合导航系统的欺骗干扰检测装置，其特征在于：所述第一处理单元具体被配置为：

利用所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值，作为所述检测统计量，构造广义似然比公式，得到零速检测公式：其中， 分别为所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值的测量方差，Zn为第n次检测的检测统计量集合，所述检测统计量为所述测量惯性单元IMU在检测时间段内各时刻的测量值，具体包括第k测量时刻的三轴加速度 和角速度N为检测窗口内数据个数，g为重力常量， 为检测窗口长度内的加速度平均值， 为加速度平均值 的幅值，T表示矩阵转置，Ωn为检测窗口长度， ， 表示

的二范数运算， 表示 的二范数运算；

其中，当所述检测统计量T(zn)小于设定的门限γ时，确定所述载体处于所述静止状态；

所述第二处理单元具体被配置为：

所述组合导航系统采用所述卡尔曼滤波模型，第k个测量时刻的状态向量为：T

xk＝[δpk，δvk，δφk，δωk，δfk，bclk]其中， 表示位置误差向量，由经度、纬度和高度组成；δvk＝[δvE，k，δvN，k，δvC，k]表示地球参考速度误差向量；δφk＝[δpk，δrk，δAk]为姿态误差向量；δωk＝[δωR，k，δωF，k，δωU，k]为陀螺仪零偏误差向量；

δfk＝[δfR，k，δfF，k，δfU，k]为加速度计零偏误差向量；bclk＝[δbr，δdr]为接收机钟差和钟漂；

其中，下标E,N,C分别表示当地地理坐标系中的东、北、天分量，下标R，F，U分别表示载体坐标系中的右、前、上分量；第k个测量时刻的测量向量zk表示如下：其中， 分别为k时刻第i颗卫星的观测伪距和观测伪距率，i＝1,2,L,M，M为观测卫星数量；则系统状态方程和测量方程为：其中，xk和xk‑1分别表示k和k‑1时刻的状态向量，wk‑1表示k‑1时刻的状态噪声向量，Φk/k‑1为k‑1时刻的状态转移矩阵，由所述惯性导航系统INS的误差方程和所述全球卫星导航系统GNSS接收机的钟差模型得出，连续到离散过程采用一阶近似；Hk表示k时刻测量矩阵，vk为k时刻观测噪声向量，wk为k时刻状态噪声向量，wk和vk皆假设为零均值的高斯白噪声，协方差分别用Qk和Rk表示；

在卡尔曼滤波的基础上，利用当前时刻的测量向量Zk和状态预测向量 构建新息向量δZk为：

7.根据权利要求6所述的一种组合导航系统的欺骗干扰检测装置，其特征在于，所述第三处理单元具体被配置为：

当检测到所述载体处于所述静止状态时，所述载体的车体坐标系b系下三个方向的速b

度为零，v＝0，在该k时刻由所述全球卫星导航系统GNSS接收机解算得到的所述载体的速度测量值为 对所述全球卫星导航系统GNSS接收机进行欺骗干扰检测，利用所述全球卫星导航系统GNSS接收机解算得到的所述载体的速度测量值作为所述静止状态下的欺骗检测统计量，构造广义似然比公式，得到所述静止状态下的欺骗干扰检测公式：其中，zn为第n次检测的检测统计量集合，N为检测窗口内数据个数，Ωn为检测窗口长度， 为GNSS观测方差， 表示的二范数运算；

当所述静止状态下的欺骗检测统计量T(zn)小于设定的门限γ时，确定所述全球卫星导航系统GNSS接收机未受到欺骗干扰，获取未受欺骗干扰时段的新息序列δzk作为所述参考噪声，用于所述运动状态下的欺骗干扰检测，所述未受欺骗干扰时段的新息序列的均值为零，统计不同卫星新息序列的方差为：其中， 为第j颗卫星新息的方差，M为统计新息个数，ΩM为未受欺骗干扰的窗口长度，δzk，j为k时刻第j颗卫星的新息。

8.根据权利要求7所述的一种组合导航系统的欺骗干扰检测装置，其特征在于，所述第四处理单元具体被配置为：

在所述静止状态下得确定所述参考噪声方差，并且当所述载体处于所述运动状态时，将欺骗速率作为满足非零均值高斯分布的随机变量，设所述运动状态下k时刻第j颗可见卫星的伪距率新息为 当所述全球卫星导航系统GNSS接收机未受到所述欺骗干扰时，伪距率新息方差应为所述静止状态下统计的对应参考噪声方差，当所述全球卫星导航系统GNSS接收机受到所述欺骗干扰时，欺骗速率的抖动使得新息方差大于参考噪声方差；令H0为无欺骗状态，H1为欺骗状态，则所述运动状态下的欺骗干扰检测为二元假设模型：

其中，检测序列长度为M，w[k]为已知参考噪声方差 的零均值高斯噪声，s[k]为由欺骗信号引入的方差 未知的零均值高斯分布信号；则广义似然比为：采用MLE方法估计 通过取对数得到：对 求导得到：

令导数为0，得到欺骗信号方差的MLE估计为：若 则相应的MLE为 与参数约束保持一致，则 的MLE为令 其中上标+表示

为正时，其为MLE；

当 时，判定为H1，将估计值 代入对数似然比可得：

‑1

令 g(z)＝z‑ln  z‑1为单调递增函数，对于z＞1，逆g 存在，若或 判定为H1，欺骗检测统计量

为对欺骗信号方差的MLE估计值

2

在H0条件下，新息序列满足零均值方差为σj的高斯分布，当样本数M足够时，新息方差渐进满足高斯分布：

其中，

则 的概率密度函数为：

根据 的概率密度函数和虚警概率计算门限，门限γ为：‑1

其中，Q 为逆高斯右尾概率函数，PFA为虚警概率，当统计量 大于设定的门限γ时判定所述载体的所述全球卫星导航系统GNSS接收机受到欺骗干扰。

9.一种电子设备，其特征在于，所述电子设备包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时，实现权利要求1‑4任一项所述的一种组合导航系统的欺骗干扰检测方法。

10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时，实现权利要求1‑4任一项所述的一种组合导航系统的欺骗干扰检测方法。

说明书

一种组合导航系统的欺骗干扰检测方法与装置

技术领域

[0001] 本发明属于卫星导航技术领域，尤其涉及一种组合导航系统的欺骗干扰检测方法与装置。

背景技术

[0002] GNSS欺骗干扰逐渐成为GNSS民用用户不容忽视的威胁，其中生成式欺骗干扰由于时延可根据需要动态调整使得接收机较难识别和消除，成为目前欺骗式干扰检测技术的研
究重点。在欺骗式干扰检测技术中，通过外部设备辅助GNSS进行欺骗检测的技术，由于引入
外部设备，使得导航信息冗余而具有良好的检测性能，常用的辅助设备是不受无线电信号
干扰的INS，其与GNSS通过卡尔曼滤波等方式进行融合并结合检测理论实现对欺骗干扰的
检测。但是传统的基于GNSS/INS组合实施欺骗检测的理论认为欺骗信号速率为一未知常
数，实际情况中由于环境噪声和各种误差等因素的影响，对目标接收机施加欺骗信号的速
率会存在抖动，将其当作常数是不合理的。

发明内容

[0003] 为解决现有技术中存在的问题，本发明提出GNSS(Global Navigation Satellite System,全球卫星导航系统)与INS(Inertial Navigation System,惯性导航系统)组合导
航系统模式下，利用测量信号进行GNSS接收机欺骗干扰检测的方案，特别适用于非连续的
缓变式欺骗干扰检测。

[0004] 本发明将存在抖动的欺骗速率建模为服从非零均值高斯分布的随机变量，首先根据INS输出的加速度和角速度测量值进行零速检测，当检测到载体静止时，此时判断GNSS信
号是否被欺骗是很容易的，根据GNSS接收机解算的用户速度是否发生变化即可判断。如果
静止时检测到欺骗直接告警，当静止状态下没有检测出欺骗时，收集此时的新息并统计方
差。由于未受欺骗时的新息方差只包括INS传播误差和观测噪声，可将此时的新息视作参考
噪声，统计参考噪声的方差作为无欺骗状态下的先验信息。当载体运动时统计一段时间内
新息的方差并与先验参考噪声方差进行比较，将欺骗速率当作随机变量处理将导致欺骗存
在时的新息方差增大，通过比较运动状态下统计的新息序列方差和静止状态下参考噪声方
差即可检测欺骗。若检测出欺骗，将欺骗信号从组合导航系统中剔除；若未检测出欺骗，则
利用GNSS信号对惯导进行校正后继续执行下一周期的流程。

[0005] 本发明第一方面公开了一种组合导航系统的欺骗干扰检测方法。所述组合导航系统包括全球卫星导航系统GNSS接收机和惯性导航系统INS，所述惯性导航系统INS包括测量
惯性单元IMU；所述方法包括：步骤S1、获取所述测量惯性单元IMU输出的三维加速度和角速
度测量值，作为检测统计量，并基于所述检测统计量构造零速检测统计量，所述零速检测统
计量用于确定载体是否处于静止状态；步骤S2、构造滤波模型，所述滤波模型用于估计当前
系统状态，并构造欺骗检测统计量，其中，所述滤波模型为卡尔曼滤波模型；步骤S3、响应于所述载体处于所述静止状态，根据所述当前系统状态构造所述静止状态下的欺骗检测统计
量，所述静止状态下的欺骗检测统计量用于确定所述载体在所述静止状态下是否受到欺骗
干扰，并且响应于所述载体在所述静止状态下未受到所述欺骗干扰，进一步确定参考噪声；
步骤S4、计算所述参考噪声的方差，利用所述参考噪声的方差构造运动状态下的欺骗干扰
统计量，所述运动状态下的欺骗干扰统计量用于判断所述载体在所述运动状态下是否受到
欺骗干扰；其中，所述载体为车辆。

[0006] 根据本发明第一方面的方法，在所述步骤S1中：

[0007] 利用所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值，作为所述检测统计量，构造广义似然比公式，得到零速检测公式：

[0008]

[0009] 其中， 分别为所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值的测量方差，Zn为第n次检测的检测统计量集合，所述检测统计量为所述测量惯性
单元IMU在检测时间段内各时刻的测量值，具体包括第k测量时刻的三轴加速度 和角速
度 N为检测窗口内数据个数，g为重力常量， 为检测窗口长度内的
加速度平均值， 为加速度平均值 的幅值，T表示矩阵转置，Ωn为检测
窗口长度， ， 表示
的二范数运算， 表示 的二范数运算；

[0010] 其中，当所述检测统计量T(zn)小于设定的门限γ时，确定所述载体处于所述静止状态。

[0011] 根据本发明第一方面的方法，在所述步骤S2中：

[0012] 所述组合导航系统采用所述卡尔曼滤波模型，第k个测量时刻的状态向量为：

[0013] xk＝[δpk，δvk，δφk，δωk，δfk，bclk]T

[0014] 其中， 表示位置误差向量，由经度、纬度和高度组成；δvk＝[δvE，k，δvN，k，δvC，k]表示地球参考速度误差向量；δφk＝[δpk，δrk，δAk]为姿态误差向量；
δωk＝[δωR，k，δωF，k，δωU，k]为陀螺仪零偏误差向量；δfk＝[δfR，k，δfF，k，δfU，k]为加速度计零偏误差向量；bclk＝[δbr，δdr]为接收机钟差和钟漂；其中，下标E，N，C分别表示当

[0015] 地地理坐标系中的东、北、天分量，下标R，F，U分别表示载体坐标系中的右、前、上分量；第k个测量时刻的测量向量zk表示如下：

[0016]

[0017] 其中， k分别为k时刻第i颗卫星的观测伪距和观测伪距率，i＝1，2，L，M，M为观测卫星数量；则系统状态方程和测量方程为：

[0018]

[0019] 其中，xk和xk‑1分别表示k和k‑1时刻的状态向量，wk‑1表示k‑1时刻的状态噪声向量，Φk/k‑1为k‑1时刻的状态转移矩阵，由所述惯性导航系统INS的误差方程和所述全球卫星导航系统GNSS接收机的钟差模型得出，连续到离散过程采用一阶近似；Hk表示k时刻测量矩阵，vk为k时刻观测噪声向量，wk为k时刻状态噪声向量，wk和vk皆假设为零均值的高斯白噪
声，协方差分别用Qk和Rk表示；

[0020] 在卡尔曼滤波的基础上，利用当前时刻的测量向量zk和状态预测向量 构建新息向量δzk为： 。

[0021] 根据本发明第一方面的方法，在所述步骤S3中：

[0022] 当检测到所述载体处于所述静止状态时，所述载体的车体坐标系b系下三个方向b
的速度为零，v ＝0，在该k时刻由所述全球卫星导航系统GNSS接收机解算得到的所述载体
的速度测量值为 对所述全球卫星导航系统GNSS接收机进行欺骗干扰检测，利用所述
全球卫星导航系统GNSS接收机解算得到的所述载体的速度测量值作为所述静止状态下的
欺骗检测统计量，构造广义似然比公式，得到所述静止状态下的欺骗干扰检测公式：

[0023]

[0024] 其中，zn为第n次检测的检测统计量集合，N为检测窗口内数据个数，Ωn为检测窗口长度， 为GNSS观测方差， 表示
的二范数运算；

[0025] 当所述静止状态下的欺骗检测统计量T(zn)小于设定的门限γ时，确定所述全球卫星导航系统GNSS接收机未受到欺骗干扰，获取未受欺骗干扰时段的新息序列δzk作为所述
参考噪声，用于所述运动状态下的欺骗干扰检测，所述未受欺骗干扰时段的新息序列的均
值为零，统计不同卫星新息序列的方差为：

[0026]

[0027] 其中， 为第j颗卫星新息的方差，M为统计新息个数，ΩM为未受欺骗干扰的窗口长度，δzk，j为k时刻第j颗卫星的新息。

[0028] 根据本发明第一方面的方法，在所述步骤S4中：

[0029] 在所述静止状态下得确定所述参考噪声方差，并且当所述载体处于所述运动状态时，将欺骗速率作为满足非零均值高斯分布的随机变量，设所述运动状态下k时刻第j颗可
见卫星的伪距率新息为 当所述全球卫星导航系统GNSS接收机未受到所述欺骗干扰
时，伪距率新息方差应为所述静止状态下统计的对应参考噪声方差，当所述全球卫星导航
系统GNSS接收机受到所述欺骗干扰时，欺骗速率的抖动使得新息方差大于参考噪声方差；
令 H0为无欺骗状态，H1为欺骗状态，则所述运动状态下的欺骗干扰检测为二
元假设模型：

[0030]

[0031] 其中，检测序列长度为M，w[k]为已知参考噪声方差 的零均值高斯噪声，s[k]为由欺骗信号引入的方差 未知的零均值高斯分布信号；则广义似然比为：

[0032]

[0033]

[0034] 采用MLE方法估计 通过取对数得到：

[0035]

[0036] 对 求导得到：

[0037]

[0038] 令导数为0，得到欺骗信号方差的MLE估计为：

[0039]

[0040] 若 ，则相应的MLE为 ，与参数约束保持一致，则 的MLE为令 其中上标+表示
为正时，其为MLE；

[0041] 当 时，判定为H1，将估计值 代入对数似然比可得：

[0042]

[0043] 令 g(z)＝z‑ln z‑1为单调递增函数，对于z＞1，逆g‑1存在，若或 判定为H1，欺骗检测统计量为
对欺骗信号方差的MLE估计值

[0044] 在H0条件下，新息序列满足零均值方差为σj2的高斯分布，当样本数M足够时，新息方差渐进满足高斯分布：

[0045]

[0046] 其 中 ，则 的概率密度函数为：

[0047]

[0048] 根据 的概率密度函数和虚警概率计算门限，门限γ为：

[0049]‑1

[0050] 其中，Q 为逆高斯右尾概率函数，PFA为虚警概率，当统计量 大于设定的门限γ时判定所述载体的所述全球卫星导航系统GNSS接收机受到欺骗干扰。

[0051] 本发明第二方面公开了一种组合导航系统的欺骗干扰检测装置。所述组合导航系统包括全球卫星导航系统GNSS接收机和惯性导航系统INS，所述惯性导航系统INS包括测量
惯性单元IMU；所述装置包括：第一处理单元，被配置为：获取所述测量惯性单元IMU输出的
三维加速度和角速度测量值，作为检测统计量，并基于所述检测统计量构造零速检测统计
量，所述零速检测统计量用于确定载体是否处于静止状态；第二处理单元，被配置为：构造
滤波模型，所述滤波模型用于估计当前系统状态，并构造欺骗检测统计量，其中，所述滤波
模型为卡尔曼滤波模型；第三处理单元，被配置为：响应于所述载体处于所述静止状态，根
据所述当前系统状态构造所述静止状态下的欺骗检测统计量，所述静止状态下的欺骗检测
统计量用于确定所述载体在所述静止状态下是否受到欺骗干扰，并且响应于所述载体在所
述静止状态下未受到所述欺骗干扰，进一步确定参考噪声；第四处理单元，被配置为：计算
所述参考噪声的方差，利用所述参考噪声的方差构造运动状态下的欺骗干扰统计量，所述
运动状态下的欺骗干扰统计量用于判断所述载体在所述运动状态下是否受到欺骗干扰；其
中，所述载体为车辆。

[0052] 根据本发明第二方面的装置，所述第一处理单元具体被配置为：

[0053] 利用所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值，作为所述检测统计量，构造广义似然比公式，得到零速检测公式：

[0054]

[0055] 其中， 分别为所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值的测量方差，zn为第n次检测的检测统计量集合，所述检测统计量为所述测量惯性单
元IMU在检测时间段内各时刻的测量值，具体包括第k测量时刻的三轴加速度 和角速度
N为检测窗口内数据个数，g为重力常量， 为检测窗口长度内的加速
度平均值， 为加速度平均值 的幅值，T表示矩阵转置，Ωn为检测窗口长
度， 表示
的二范数运算， 表示 的二范数运算；

[0056] 其中，当所述检测统计量T(zn)小于设定的门限γ时，确定所述载体处于所述静止状态。

[0057] 根据本发明第二方面的装置，所述第二处理单元具体被配置为：

[0058] 所述组合导航系统采用所述卡尔曼滤波模型，第k个测量时刻的状态向量为：

[0059] xk＝[δpk，δvk，δφk，δωk，δfk，bclk]T

[0060] 其中， 表示位置误差向量，由经度、纬度和高度组成；δvk＝[δvE，k，δvN，k，δvC，k]表示地球参考速度误差向量；δφk＝[δpk，δrk，δAk]为姿态误差向量；
δωk＝[δωR，k，δωF，k，δωU，k]为陀螺仪零偏误差向量；δfk＝[δfR，k，δfF，k，δfU，k]为加速度计零偏误差向量；bclk＝[δbr，δdr]为接收机钟差和钟漂；其中，下标E，N，C分别表示当地地理坐标系中的东、北、天分量，下标R，F，U分别表示载体坐标系中的右、前、上分量；第k个测量时刻的测量向量zk表示如下：

[0061]

[0062] 其中， k分别为k时刻第i颗卫星的观测伪距和观测伪距率，i＝1，2，L，M，M为观测卫星数量；则系统状态方程和测量方程为：

[0063]

[0064] 其中，xk和xk‑1分别表示k和k‑1时刻的状态向量，wk‑1表示k‑1时刻的状态噪声向量，Φk/k‑1为k‑1时刻的状态转移矩阵，由所述惯性导航系统INS的误差方程和所述全球卫星导航系统GNSS接收机的钟差模型得出，连续到离散过程采用一阶近似；Hk表示k时刻测量矩阵，vk为k时刻观测噪声向量，wk为k时刻状态噪声向量，wk和vk皆假设为零均值的高斯白噪
声，协方差分别用Qk和Rk表示；

[0065] 在卡尔曼滤波的基础上，利用当前时刻的测量向量zk和状态预测向量 构建新息向量δzk为：

[0066] 根据本发明第二方面的装置，所述第三处理单元具体被配置为：

[0067] 当检测到所述载体处于所述静止状态时，所述载体的车体坐标系b系下三个方向b
的速度为零，v ＝0，在该k时刻由所述全球卫星导航系统GNSS接收机解算得到的所述载体
的速度测量值为 对所述全球卫星导航系统GNSS接收机进行欺骗干扰检测，利用所述
全球卫星导航系统GNSS接收机解算得到的所述载体的速度测量值作为所述静止状态下的
欺骗检测统计量，构造广义似然比公式，得到所述静止状态下的欺骗干扰检测公式：

[0068]

[0069] 其中，zn为第n次检测的检测统计量集合，N为检测窗口内数据个数，Ωn为检测窗口长度， 为GNSS观测方差， 表示
的二范数运算；

[0070] 当所述静止状态下的欺骗检测统计量T(zn)小于设定的门限γ时，确定所述全球卫星导航系统GNSS接收机未受到欺骗干扰，获取未受欺骗干扰时段的新息序列δzk作为所述
参考噪声，用于所述运动状态下的欺骗干扰检测，所述未受欺骗干扰时段的新息序列的均
值为零，统计不同卫星新息序列的方差为：

[0071]

[0072] 其中， 为第j颗卫星新息的方差，M为统计新息个数，ΩM为未受欺骗干扰的窗口长度，δzk，j为k时刻第j颗卫星的新息。

[0073] 根据本发明第二方面的装置，所述第四处理单元具体被配置为：

[0074] 在所述静止状态下得确定所述参考噪声方差，并且当所述载体处于所述运动状态时，将欺骗速率作为满足非零均值高斯分布的随机变量，设所述运动状态下k时刻第j颗可
见卫星的伪距率新息为 当所述全球卫星导航系统GNSS接收机未受到所述欺骗干扰
时，伪距率新息方差应为所述静止状态下统计的对应参考噪声方差，当所述全球卫星导航
系统GNSS接收机受到所述欺骗干扰时，欺骗速率的抖动使得新息方差大于参考噪声方差；
令 H0为无欺骗状态，H1为欺骗状态，则所述运动状态下的欺骗干扰检测为二
元假设模型：

[0075]

[0076] 其中，检测序列长度为M，w[k]为已知参考噪声方差 的零均值高斯噪声，s[k]为由欺骗信号引入的方差 未知的零均值高斯分布信号；则广义似然比为：

[0077]

[0078]

[0079] 采用MLE方法估计 通过取对数得到：

[0080]

[0081] 对 求导得到：

[0082]

[0083] 令导数为0，得到欺骗信号方差的MLE估计为：

[0084]

[0085] 若 则相应的MLE为 与参数约束保持一致，则 的MLE为令 其中上标+表示
为正时，其为MLE；

[0086] 当 时，判定为H1，将估计值 代入对数似然比可得：

[0087]‑1

[0088] 令 g(z)＝z‑ln z‑1为单调递增函数，对于z＞1，逆g 存在，若或 判定为H1，欺骗检测统计量
为对欺骗信号方差的MLE估计值

[0089] 在H0条件下，新息序列满足零均值方差为σj2的高斯分布，当样本数M足够时，新息方差渐进满足高斯分布：

[0090]

[0091] 其中，则 的概率密度函数为：

[0092]

[0093] 根据 的概率密度函数和虚警概率计算门限，门限γ为：

[0094]‑1

[0095] 其中，Q 为逆高斯右尾概率函数，PFA为虚警概率，当统计量 大于设定的门限γ时判定所述载体的所述全球卫星导航系统GNSS接收机受到欺骗干扰。

[0096] 本发明第三方面公开了一种电子设备。所述电子设备包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时，实现本公开第一方面所述
的一种组合导航系统的欺骗干扰检测方法中的步骤。

[0097] 本发明第四方面公开了一种计算机可读存储介质。所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时，实现本公开第一方面所述的一种组合
导航系统的欺骗干扰检测方法中的步骤。

[0098] 综上，本发明提出的技术方案本发明将抖动的欺骗速率建模为服从非零均值高斯分布的随机变量，通过比较静止无欺骗状态下统计的参考噪声方差和运动状态下的新息序
列方差，构建新息序列方差的检测统计量，进行欺骗式干扰检测。由于设计的检测方案主要
检测方差的变化，所以其对欺骗速率均值不敏感，对检测速率缓变的欺骗干扰具有很好的
检测性能。本发明带来的明显技术效果主要如下：(1)提出了针对抖动欺骗速率的新息序列
方差欺骗干扰检测方法，相比传统的利用新息均值变化进行欺骗检测的SPRT(Sequential 
Probability Ratio Test，序贯概率比检验)方法，该方法对缓变式欺骗干扰检测的灵敏度
较高，能够及时检出欺骗的产生和消失，特别适用于非连续的缓变式欺骗干扰检测，可提高
GNSS信号的利用率。(2)在多星受到欺骗时，该方法能够较好地识别出受欺骗卫星，同时未
被欺骗的卫星信号对应的检测统计量不会受到被欺骗卫星信号的影响。且欺骗检测率在欺
骗速率方差较小时已经较高，受欺骗速率均值变化影响较小，对检测速率缓变的欺骗干扰
具有较好的检测率性能。

附图说明

[0099] 为了更清楚地说明本发明具体实施方式或现有技术中的技术方案下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的
附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前
提下，还可以根据这些附图获得其他的附图。

[0100] 图1为根据本发明实施例的一种组合导航系统的欺骗干扰检测方法的流程示意图；

[0101] 图2a‑2b为本发明与传统SPRT算法对均值0.2m/s，方差0.1(m/s)2欺骗速率的检测灵敏度对比的示意图；

[0102] 图3a‑3b为本发明与传统SPRT算法对均值0.4m/s，方差0.1(m/s)2欺骗速率的检测灵敏度对比的示意图；

[0103] 图4a‑4b为本发明与传统SPRT算法对均值0.4m/s，方差0.05(m/s)2欺骗速率的检测灵敏度对比的示意图；

[0104] 图5为根据本发明实施例的对多星进行欺骗检测的结果图；

[0105] 图6为根据本发明实施例的欺骗检测率与欺骗速率方差和均值的关系图；

[0106] 图7为根据本发明实施例的一种电子设备的结构图。

具体实施方式

[0107] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例只
是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人
员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0108] 本发明第一方面公开了一种组合导航系统的欺骗干扰检测方法。所述组合导航系统包括全球卫星导航系统GNSS接收机和惯性导航系统INS，所述惯性导航系统INS包括测量
惯性单元IMU；所述方法包括：步骤S1、获取所述测量惯性单元IMU输出的三维加速度和角速
度测量值，作为检测统计量，并基于所述检测统计量构造零速检测统计量，所述零速检测统
计量用于确定载体是否处于静止状态；步骤S2、构造滤波模型，所述滤波模型用于估计当前
系统状态，并构造欺骗检测统计量，其中，所述滤波模型为卡尔曼滤波模型；步骤S3、响应于所述载体处于所述静止状态，根据所述当前系统状态构造所述静止状态下的欺骗检测统计
量，所述静止状态下的欺骗检测统计量用于确定所述载体在所述静止状态下是否受到欺骗
干扰，并且响应于所述载体在所述静止状态下未受到所述欺骗干扰，进一步确定参考噪声；
步骤S4、计算所述参考噪声的方差，利用所述参考噪声的方差构造运动状态下的欺骗干扰
统计量，所述运动状态下的欺骗干扰统计量用于判断所述载体在所述运动状态下是否受到
欺骗干扰；其中，所述载体为车辆。

[0109] 在一些实施例中，在所述步骤S1中：

[0110] 利用所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值，作为所述检测统计量，构造广义似然比公式，得到零速检测公式：

[0111]

[0112] 其中， 分别为所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值的测量方差，zn为第n次检测的检测统计量集合，所述检测统计量为所述测量惯性单
元IMU在检测时间段内各时刻的测量值，具体包括第k测量时刻的三轴加速度 和角速度
N为检测窗口内数据个数，g为重力常量， 为检测窗口长度内的
加速度平均值， 为加速度平均值 的幅值，T表示矩阵转置，Ωn为检测
窗口长度， ， 表示
的二范数运算， 表示 的二范数运算；

[0113] 其中，当所述检测统计量T(zn)小于设定的门限γ时，确定所述载体处于所述静止状态。

[0114] 在一些实施例中，在所述步骤S2中：

[0115] 所述组合导航系统采用所述卡尔曼滤波模型，第k个测量时刻的状态向量为：

[0116] xk＝[δpk，δvk，δφk，δωk，δfk，bclk]T

[0117] 其中， 表示位置误差向量，由经度、纬度和高度组成；δvk＝[δvE，k，δvN，k，δvC，k]表示地球参考速度误差向量；δφk＝[δpk，δrk，δAk]为姿态误差向量；
δωk＝[δωR，k，δωF，k，δωU，k]为陀螺仪零偏误差向量；δfk＝[δfR，k，δfF，k，δfU，k]为加速度计零偏误差向量；bclk＝[δbr，δdr]为接收机钟差和钟漂；其中，下标E，N，C分别表示当地地理坐标系中的东、北、天分量，下标R，F，U分别表示载体坐标系中的右、前、上分量；第k个测量时刻的测量向量zk表示如下：

[0118]

[0119] 其中， k分别为k时刻第i颗卫星的观测伪距和观测伪距率，i＝1，2，L，M，M为观测卫星数量；则系统状态方程和测量方程为：

[0120]

[0121] 其中，xk和xk‑1分别表示k和k‑1时刻的状态向量，wk‑1表示k‑1时刻的状态噪声向量，Φk/k‑1为k‑1时刻的状态转移矩阵，由所述惯性导航系统INS的误差方程和所述全球卫星导航系统GNSS接收机的钟差模型得出，连续到离散过程采用一阶近似；Hk表示k时刻测量矩阵，vk为k时刻观测噪声向量，wk为k时刻状态噪声向量，wk和vk皆假设为零均值的高斯白噪
声，协方差分别用Qk和Rk表示；

[0122] 在卡尔曼滤波的基础上，利用当前时刻的测量向量zk和状态预测向量 构建新息向量δzk为：

[0123] 在一些实施例中，在所述步骤S3中：

[0124] 当检测到所述载体处于所述静止状态时，所述载体的车体坐标系b系下三个方向b
的速度为零，v ＝0，在该k时刻由所述全球卫星导航系统GNSS接收机解算得到的所述载体
的速度测量值为 ，对所述全球卫星导航系统GNSS接收机进行欺骗干扰检测，利用所
述全球卫星导航系统GNSS接收机解算得到的所述载体的速度测量值作为所述静止状态下
的欺骗检测统计量，构造广义似然比公式，得到所述静止状态下的欺骗干扰检测公式：

[0125]

[0126] 其中，zn为第n次检测的检测统计量集合，N为检测窗口内数据个数，Ωn为检测窗口长度， 为GNSS观测方差， 表示
的二范数运算；

[0127] 当所述静止状态下的欺骗检测统计量T(zn)小于设定的门限γ时，确定所述全球卫星导航系统GNSS接收机未受到欺骗干扰，获取未受欺骗干扰时段的新息序列δzk作为所述
参考噪声，用于所述运动状态下的欺骗干扰检测，所述未受欺骗干扰时段的新息序列的均
值为零，统计不同卫星新息序列的方差为：

[0128]

[0129] 其中， 为第j颗卫星新息的方差，M为统计新息个数，ΩM为未受欺骗干扰的窗口长度，δzk，j为k时刻第j颗卫星的新息。

[0130] 在一些实施例中，在所述步骤S4中：

[0131] 在所述静止状态下得确定所述参考噪声方差，并且当所述载体处于所述运动状态时，将欺骗速率作为满足非零均值高斯分布的随机变量，设所述运动状态下k时刻第j颗可
见卫星的伪距率新息为 当所述全球卫星导航系统GNSS接收机未受到所述欺骗干扰时，
伪距率新息方差应为所述静止状态下统计的对应参考噪声方差，当所述全球卫星导航系统
GNSS接收机受到所述欺骗干扰时，欺骗速率的抖动使得新息方差大于参考噪声方差；令
H0为无欺骗状态，H1为欺骗状态，则所述运动状态下的欺骗干扰检测为二元
假设模型：

[0132]

[0133] 其中，检测序列长度为M，w[k]为已知参考噪声方差 的零均值高斯噪声，s[k]为由欺骗信号引入的方差 未知的零均值高斯分布信号；则广义似然比为：

[0134]

[0135]

[0136] 采用MLE方法估计 通过取对数得到：

[0137]

[0138] 对 求导得到：

[0139]

[0140] 令导数为0，得到欺骗信号方差的MLE估计为：

[0141]

[0142] 若 则相应的MLE为 与参数约束保持一致，则 的MLE为令 其中上标+表示
为正时，其为MLE；

[0143] 当 时，判定为H1，将估计值 代入对数似然比可得：

[0144]‑1

[0145] 令 g(z)＝z‑ln z‑1为单调递增函数，对于z＞1，逆g 存在，若或 判定为H1，欺骗检测统计量为对
欺骗信号方差的MLE估计值

[0146] 在H0条件下，新息序列满足零均值方差为σj2的高斯分布，当样本数M足够时，新息方差渐进满足高斯分布：

[0147]

[0148] 其中，则 的概率密度函数为：

[0149]

[0150] 根据 的概率密度函数和虚警概率计算门限，门限γ为：

[0151]

[0152] 其中，Q‑1为逆高斯右尾概率函数，PFA为虚警概率，当统计量 大于设定的门限γ时判定所述载体的所述全球卫星导航系统GNSS接收机受到欺骗干扰。

[0153] 结合图1(包括全球卫星导航系统GNSS和惯性导航系统INS的组合导航系统的欺骗干扰检测方法的流程图)对所述方法进行详细说明。

[0154] 第一步，构造零速检测统计量：

[0155] 作为一个更优的实施例，该步骤利用IMU输出的三维加速度、角速度测量值作为检测统计量，构造广义似然比公式，得到如下零速检测公式：

[0156]

[0157] 其中， 分别为所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值的测量方差，zn为第n次检测的检测统计量集合，所述检测统计量为所述测量惯性单
元IMU在检测时间段内各时刻的测量值，具体包括第k测量时刻的三轴加速度 和角速度
N为检测窗口内数据个数，g为重力常量， 为检测窗口长度内的加
速度平均值， 为加速度平均值 的幅值，T表示矩阵转置，Ωn为检测窗
口长度， 表示
的二范数运算， 表示 的二范数运算；

[0158] 当所述检测统计量T(zn)小于设定的门限γ时认为载体处于静止状态。

[0159] 第二步，构造滤波模型：

[0160] 滤波模型采用常用的卡尔曼滤波：卡尔曼滤波模型包括状态方程和测量方程；状态方程由INS的误差传播方程和GNSS接收机钟差模型得出，并对连续过程采用离散一阶近
似；测量方程由GNSS观测模型得出。在任意一个测量时刻i，卡尔曼滤波系统的状态量为xi，对该时刻状态量的估计值由预测值和观测值经卡尔曼增益加权后得出。

[0161] 作为一个更优的实施方式，该步骤中的具体包括：

[0162] 第k个测量时刻的状态向量为：

[0163] xk＝[δpk，δvk，δφk，δωk，δfk，bclk]T

[0164] 其中， 表示位置误差向量，由经度、纬度和高度组成；δvk＝[δvE，k，δvN，k，δvC，k]表示地球参考速度误差向量；δφk＝[δPk，δrk，δAk]为姿态误差向量；
δωk＝[δωR，k，δωF，k，δωU，k]为陀螺仪零偏误差向量；δfk＝[δfR，k，δfF，k，δfU，k]为加速度计零偏误差向量；bclk＝[δbr，δdr]为接收机钟差和钟漂；其中，下标E，N，C分别表示当地地理坐标系中的东、北、天分量，下标R，F，U分别表示载体坐标系中的右、前、上分量；

[0165] 第k个测量时刻的测量向量zk表示如下：

[0166]

[0167] 其中， (i＝1，2，L，M，M为观测卫星数量)分别为k时刻第i颗卫星的观测伪距和观测伪距率。

[0168] 系统的状态方程和测量方程可写为：

[0169]

[0170] 其中，xk和xk‑1分别表示k和k‑1时刻的状态向量，wk‑1表示k‑1时刻的状态噪声向量，Φk/k‑1为k‑1时刻的状态转移矩阵，由所述惯性导航系统INS的误差方程和所述全球卫星导航系统GNSS接收机的钟差模型得出，连续到离散过程采用一阶近似；Hk表示k时刻测量矩阵，vk为k时刻观测噪声向量，wk为k时刻状态噪声向量，wk和vk皆假设为零均值的高斯白噪
声，协方差分别用Qk和Rk表示；

[0171] 在卡尔曼滤波的基础上，利用当前时刻的测量向量zk和状态预测向量 可构建新息向量δzk为：

[0172]

[0173] 第三步，构造静止状态欺骗检测统计量：

[0174] 当检测到载体静止时，即载体的车体坐标系b系下三个方向的速度为零，即vb＝0，而在该k时刻由GNSS接收机解算得到的载体速度测量值为 ，此时可进行GNSS接收机欺
骗干扰检测，利用接收机解算后输出的载体速度作为检测统计量，构造广义似然比公式，得
到如下静止状态欺骗干扰检测公式：

[0175]

[0176] 其中，zn为第n次检测的检测统计量集合，N为检测窗口内数据个数，Ωn为检测窗口长度， 为GNSS观测方差， 表示
的二范数运算；

[0177] 当所述检测统计量T(zn)小于设定的门限γ时认为GNSS接收机未受到欺骗干扰，获取未受欺骗干扰时段的新息序列δzk作为参考噪声用于运动状态欺骗干扰检测，由于未受
欺骗条件下新息的均值为零，统计不同卫星新息序列的方差为：

[0178]

[0179] 其中， 为第j颗卫星新息的方差，M为统计新息个数，ΩM为未受欺骗干扰的窗口长度，δzk，j为k时刻第j颗卫星的新息。

[0180] 第四步，构造运动状态欺骗检测统计量：

[0181] 在静止状态下得到参考噪声方差后，当载体处于运动状态时，由于实际欺骗干扰速率无法保持恒定，将欺骗速率视为满足非零均值高斯分布的随机变量，设运动状态下k时
刻第j颗可见卫星的伪距率新息为 当GNSS接收机未受到欺骗干扰时，伪距率新息方差
应为静止状态下统计的对应参考噪声方差，当GNSS接收机受到欺骗干扰时，由于欺骗速率
的抖动使得新息方差将大于参考噪声方差；令 H0为无欺骗状态，H1为欺骗状
态，则所述运动状态下的欺骗干扰检测为二元假设模型：

[0182]

[0183] 其中M为检测序列长度，w[k]为已知参考噪声方差 的零均值高斯噪声，s[k]为由欺骗信号引入的方差 未知的零均值高斯分布信号；

[0184] 那么广义似然比为：

[0185]

[0186]

[0187] 采用MLE方法估计 对(2)式取对数可得：

[0188]

[0189] 将上式对 求导得：

[0190]

[0191] 令该导数为0可得欺骗信号方差的MLE估计为：

[0192]

[0193] 但是如果 那么相应的MLE应为 这与参数的约束是一致的，因此，的MLE为：

[0194]

[0195] 令 其中上标“+”表示如果 是正的，那么它是MLE；

[0196] 当 那么判H1，将估计值 代入对数似然比可得：

[0197]

[0198] 令 g(z)＝z‑ln z‑1为单调递增函数，对于z＞1，逆g‑1存在，若或 判定为H1，欺骗检测统计量为
对欺骗信号方差的MLE估计值

[0199] 在H0条件下，新息序列满足零均值方差为σj2的高斯分布，当样本数M足够时，新息方差渐进满足高斯分布：

[0200]

[0201] 其中，则 的概率密度函数为：

[0202]

[0203] 在实际应用时，虚警概率往往很小，所以在根据 的概率密度函数和虚警概率计算门限时，将 完全视为高斯分布对门限的计算没有影响；最后门限γ计算可得：

[0204]

[0205] 其中Q‑1为逆高斯右尾概率函数，PFA为虚警概率，当所述统计量 大于设定的门限γ时认为载体GNSS接收机受到欺骗干扰。

[0206] 针对本发明设计的欺骗干扰方法，进行三组仿真实验进行验证，分别为：

[0207] 实验1：在单颗卫星受到均值0.2m/s、方差0.1(m/s)2；均值0.4m/s、方差0.1(m/s)2；2
均值0.4m/s、方差0.05(m/s) 欺骗速率的缓变式欺骗时比较新息序列SPRT检测算法和本发
明的检测时间以及欺骗干扰消失后两种方法回到无欺骗状态下所需的重置时间，验证本发
明对缓变式欺骗干扰的检测灵敏度。同时观察了欺骗速率均值和方差设置对两种算法的影
响。

[0208] 实验2：在多颗卫星受到均值0.4m/s、方差0.1(m/s)2欺骗速率的缓变式欺骗时，验证本发明对多星欺骗的检测性能。

[0209] 实验3：改变欺骗速率的均值和方差，验证本发明的欺骗检测率随欺骗速率均值和方差的变化关系。

[0210] 实验相关参数设置如表1所示。

[0211] 表1

[0212]

[0213] 实验中遮蔽角设置为10°，保持可见卫星数为8颗，在实验1和3中对PRN01卫星施加欺骗，在实验2中对PRN01和PRN04卫星同时施加欺骗。

[0214] 实验1和2中缓变欺骗量的设置见表2。

[0215] 表2

[0216]

[0217] 对8颗可见卫星的伪距率参考噪声方差统计值见表3。

[0218] 表3

[0219] 2可见卫星 参考噪声方差(m/s)
PRN01 0.0117
PRN02 0.0121
PRN03 0.0137
PRN04 0.0123
PRN05 0.0132
PRN06 0.0128
PRN07 0.0123
PRN08 0.0127

[0220] 在50s开始对卫星PRN01施加均值0.2m/s、方差0.1(m/s)2；均值0.4m/s、方差0.1(m/2 2
s) ；均值0.4m/s、方差0.05(m/s)欺骗速率的缓变式欺骗，SPRT算法与本发明的欺骗检验统
计量随时间的变化曲线如图2a‑2b、图3a‑3b和图4a‑4b所示。由以上三组图可看出本发明对欺骗干扰的敏感性较强，当150s后欺骗干扰消失时，本发明的检验统计量分别在16s、16s和
15s内低于检测阈值。而对于新息序列SPRT算法而言，当欺骗干扰消失后，其检测结果依然
2
为欺骗，该状态一直持续到检测结束。此外，SPRT算法对0.2m/s均值、0.1(m/s) 方差欺骗速
2 2
率的检测时间最长，为27s；对0.4m/s均值、0.1(m/s) 和0.05(m/s)方差欺骗速率的检测时
间相同，为10s。欺骗速率的均值对SPRT算法的检测速率有一定影响。本发明对0.2m/s均值、
2 2
0.1(m/s) 方差和0.4m/s均值、0.1(m/s)方差欺骗速率的检测时间为7s，比SPRT算法快20s；
2
对0.4m/s均值、0.05(m/s) 方差欺骗速率的检测时间为11s，比SPRT算法慢1s。本发明主要
受到欺骗速率方差的影响，在实际情况中是无法对欺骗速率的方差进行控制的，欺骗速率
越大其产生的抖动也越大，因此在真实欺骗干扰环境下本发明在检测速度上具有较明显优
势。

[0221] 在50s时同时对卫星PRN01和PRN04施加均值0.4m/s、方差0.1(m/s)2欺骗速率的缓变式欺骗，本发明的检测结果如图5所示。由图可见，当对多颗卫星施加欺骗时，本发明能够较好地检测出哪些卫星受到了欺骗，同时未被欺骗的卫星信号对应的检测统计量一直处于
较低值，没有受到被欺骗卫星信号的影响。

[0222] 在50s对卫星PRN01施加均值分别为0.2m/s和0.4m/s，方差变化范围为0‑0.15(m/2
s) 欺骗速率的欺骗干扰，对所提算法的检测率进行仿直和分析，检测率P的计算公式为：

[0223]

[0224] 其中Δτ＝1s，为GNSS接收机的数据频率；Num[T＞γ，ts≤t≤ts+Δt]为欺骗干扰存在时成功检测的实验次数；Δts为欺骗干扰持续的时间。分别选择欺骗速率均值为0.2m/s和0.4m/s，分析欺骗速率方差对检测率的影响如图6所示。从图中可以看出随着欺骗速率方
差的增大，本发明的欺骗检测率逐渐提高。对于虚警概率为10‑4，在欺骗干扰方差小于0.09
2 2
(m/s) 时，检测概率小于93.07％；当方差在0.14((m/s)以上时检测概率为100％。且本方法
的检测率不受欺骗速率均值影响，特别对检测速率缓变的欺骗干扰具有很好的检测性能。

[0225] 本发明第二方面公开了一种组合导航系统的欺骗干扰检测装置。所述组合导航系统包括全球卫星导航系统GNSS接收机和惯性导航系统INS，所述惯性导航系统INS包括测量
惯性单元IMU；所述装置包括：第一处理单元，被配置为：获取所述测量惯性单元IMU输出的
三维加速度和角速度测量值，作为检测统计量，并基于所述检测统计量构造零速检测统计
量，所述零速检测统计量用于确定载体是否处于静止状态；第二处理单元，被配置为：构造
滤波模型，所述滤波模型用于估计当前系统状态，并构造欺骗检测统计量，其中，所述滤波
模型为卡尔曼滤波模型；第三处理单元，被配置为：响应于所述载体处于所述静止状态，根
据所述当前系统状态构造所述静止状态下的欺骗检测统计量，所述静止状态下的欺骗检测
统计量用于确定所述载体在所述静止状态下是否受到欺骗干扰，并且响应于所述载体在所
述静止状态下未受到所述欺骗干扰，进一步确定参考噪声；第四处理单元，被配置为：计算
所述参考噪声的方差，利用所述参考噪声的方差构造运动状态下的欺骗干扰统计量，所述
运动状态下的欺骗干扰统计量用于判断所述载体在所述运动状态下是否受到欺骗干扰；其
中，所述载体为车辆。

[0226] 在一些实施例中，所述第一处理单元具体被配置为：

[0227] 利用所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值，作为所述检测统计量，构造广义似然比公式，得到零速检测公式：

[0228]

[0229] 其中， 分别为所述测量惯性单元IMU输出的所述三维加速度和所述角速度测量值的测量方差，zn为第n次检测的检测统计量集合，所述检测统计量为所述测量惯性单
元IMU在检测时间段内各时刻的测量值，具体包括第k测量时刻的三轴加速度 和角速度
N为检测窗口内数据个数，g为重力常量， 为检测窗口长度内的加
速度平均值， 为加速度平均值 的幅值，T表示矩阵转置，Ωn为检测窗口
长度， ， 表示
的二范数运算， 表示 的二范数运算；

[0230] 其中，当所述检测统计量T(zn)小于设定的门限γ时，确定所述载体处于所述静止状态。

[0231] 在一些实施例中，述第二处理单元具体被配置为：

[0232] 所述组合导航系统采用所述卡尔曼滤波模型，第k个测量时刻的状态向量为：

[0233] xk＝[δpk，δvk，δφk，δωk，δfk，bclk]T

[0234] 其中， 表示位置误差向量，由经度、纬度和高度组成；δvx＝[δvE，k，δvN，k，δvC，k]表示地球参考速度误差向量；δφk＝[δpk，δrk，δAk]为姿态误差向量；δωk＝[δωR，k，δωF，k，δωU，k]为陀螺仪零偏误差向量；δfk＝[δfR，k，δfF，k，δfU，k]为加速度计零偏误差向量；bclk＝[δbr，δdr]为接收机钟差和钟漂；其中，下标E，N，C分别表示当地地理坐标系中的东、北、天分量，下标R，F，U分别表示载体坐标系中的右、前、上分量；第k个测量时刻的测量向量zk表示如下：

[0235]

[0236]

[0237] 其中，xk和xk‑1分别表示k和k‑1时刻的状态向量，wk‑1表示k‑1时刻的状态噪声向量，Φk/k‑1为k‑1时刻的状态转移矩阵，由所述惯性导航系统INS的误差方程和所述全球卫星导航系统GNSS接收机的钟差模型得出，连续到离散过程采用一阶近似；Hk表示k时刻测量矩阵，vk为k时刻观测噪声向量，wk为k时刻状态噪声向量，wk和vk皆假设为零均值的高斯白噪
声，协方差分别用Qk和Rk表示；

[0238] 在卡尔曼滤波的基础上，利用当前时刻的测量向量zk和状态预测向量 构建新息向量δzk为：

[0239] 在一些实施例中，所述第三处理单元具体被配置为：

[0240] 当检测到所述载体处于所述静止状态时，所述载体的车体坐标系b系下三个方向b
的速度为零，v，＝0，在该k时刻由所述全球卫星导航系统GNSS接收机解算得到的所述载体
的速度测量值为 对所述全球卫星导航系统GNSS接收机进行欺骗干扰检测，利用所述
全球卫星导航系统GNSS接收机解算得到的所述载体的速度测量值作为所述静止状态下的
欺骗检测统计量，构造广义似然比公式，得到所述静止状态下的欺骗干扰检测公式：

[0241]

[0242] 其中，zn为第n次检测的检测统计量集合，N为检测窗口内数据个数，Ωn为检测窗口长度， 为GNSS观测方差， 表示
的二范数运算；

[0243] 当所述静止状态下的欺骗检测统计量T(zn)小于设定的门限γ时，确定所述全球卫星导航系统GNSS接收机未受到欺骗干扰，获取未受欺骗干扰时段的新息序列δzk作为所述
参考噪声，用于所述运动状态下的欺骗干扰检测，所述未受欺骗干扰时段的新息序列的均
值为零，统计不同卫星新息序列的方差为：

[0244]

[0245] 其中， 为第j颗卫星新息的方差，M为统计新息个数，ΩM为未受欺骗干扰的窗口长度，δzk，j为k时刻第j颗卫星的新息。

[0246] 根据本发明第二方面的装置，所述第四处理单元具体被配置为：

[0247] 在所述静止状态下得确定所述参考噪声方差，并且当所述载体处于所述运动状态时，将欺骗速率作为满足非零均值高斯分布的随机变量，设所述运动状态下k时刻第j颗可
见卫星的伪距率新息为 当所述全球卫星导航系统GNSS接收机未受到所述欺骗干扰
时，伪距率新息方差应为所述静止状态下统计的对应参考噪声方差，当所述全球卫星导航
系统GNSS接收机受到所述欺骗干扰时，欺骗速率的抖动使得新息方差大于参考噪声方差；
令 H0为无欺骗状态，H1为欺骗状态，则所述运动状态下的欺骗干扰检测为二
元假设模型：

[0248]

[0249] 其中，检测序列长度为M，w[k]为已知参考噪声方差 的零均值高斯噪声，s[k]为由欺骗信号引入的方差 未知的零均值高斯分布信号；则广义似然比为：

[0250]

[0251]

[0252] 采用MLE方法估计 通过取对数得到：

[0253]

[0254] 对 求导得到：

[0255]

[0256] 令导数为0，得到欺骗信号方差的MLE估计为：

[0257]

[0258] 若 则相应的MLE为 与参数约束保持一致，则 的MLE为令 其中上标+表示
为正时，其为MLE；

[0259] 当 时，判定为H1，将估计值 代入对数似然比可得：

[0260]‑1

[0261] 令 g(z)＝z‑ln z‑1为单调递增函数，对于z＞1，逆g 存在，若或 判定为H1，欺骗检测统计量
为对欺骗信号方差的MLE估计值

[0262] 在H0条件下，新息序列满足零均值方差为σj2的高斯分布，当样本数M足够时，新息方差渐进满足高斯分布：

[0263]

[0264] 其中，则 的概率密度函数为：

[0265]

[0266] 根据 的概率密度函数和虚警概率计算门限，门限γ为：

[0267]

[0268] 其中，Q‑1为逆高斯右尾概率函数，PFA为虚警概率，当统计量 大于设定的门限γ时判定所述载体的所述全球卫星导航系统GNSS接收机受到欺骗干扰。

[0269] 本发明第三方面公开了一种电子设备。所述电子设备包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时，实现本公开第一方面所述
的一种组合导航系统的欺骗干扰检测方法中的步骤。

[0270] 图7为根据本发明实施例的一种电子设备的结构图，如图7所示，电子设备包括通过系统总线连接的处理器、存储器、通信接口、显示屏和输入装置。其中，该电子设备的处理器用于提供计算和控制能力。该电子设备的存储器包括非易失性存储介质、内存储器。该非
易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作
系统和计算机程序的运行提供环境。该电子设备的通信接口用于与外部的终端进行有线或
无线方式的通信，无线方式可通过WIFI、运营商网络、近场通信(NFC)或其他技术实现。该电子设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该电子设备的输入装置可以是显
示屏上覆盖的触摸层，也可以是电子设备外壳上设置的按键、轨迹球或触控板，还可以是外
接的键盘、触控板或鼠标等。

[0271] 本领域技术人员可以理解，图7中示出的结构，仅仅是与本公开的技术方案相关的部分的结构图，并不构成对本申请方案所应用于其上的电子设备的限定，具体的电子设备
可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。

[0272] 本发明第四方面公开了一种计算机可读存储介质。所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时，实现本公开第一方面所述的一种组合
导航系统的欺骗干扰检测方法中的步骤。

[0273] 综上，本发明提出的技术方案本发明将抖动的欺骗速率建模为服从非零均值高斯分布的随机变量，通过比较静止无欺骗状态下统计的参考噪声方差和运动状态下的新息序
列方差，构建新息序列方差的检测统计量，进行欺骗式干扰检测。由于设计的检测方案主要
检测方差的变化，所以其对欺骗速率均值不敏感，对检测速率缓变的欺骗干扰具有很好的
检测性能。本发明带来的明显技术效果主要如下：(1)提出了针对抖动欺骗速率的新息序列
方差欺骗干扰检测方法，相比传统的利用新息均值变化进行欺骗检测的SPRT(Sequential 
Probability Ratio Test，序贯概率比检验)方法，该方法对缓变式欺骗干扰检测的灵敏度
较高，能够及时检出欺骗的产生和消失，特别适用于非连续的缓变式欺骗干扰检测，可提高
GNSS信号的利用率。(2)在多星受到欺骗时，该方法能够较好地识别出受欺骗卫星，同时未
被欺骗的卫星信号对应的检测统计量不会受到被欺骗卫星信号的影响。且欺骗检测率在欺
骗速率方差较小时已经较高，受欺骗速率均值变化影响较小，对检测速率缓变的欺骗干扰
具有较好的检测率性能。

[0274] 请注意，以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不
存在矛盾，都应当认为是本说明书记载的范围。以上所述实施例仅表达了本申请的几种实
施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出
的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变
形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求
为准。