CN115236701A - 一种基于联合sqm平方的gnss欺骗干扰检测方法 - Google Patents

Abstract

本发明公开了一种基于联合SQM平方的GNSS欺骗干扰检测方法，包括如下步骤：获取GNSS接收机的载噪比、相关器间距、相干积分时间，并计算跟踪环路的提前路、即时路、滞后路的相干积分值；根据相干积分值计算Delta、Ratio、ELP三种SQM指标；根据接收机的载噪比、相关器间距、相干积分时间计算SQM指标的统计特征；建立联合检测量；设定虚警率，根据联合检测量的统计分布确定判决门限；设定时间窗长度，在检测时间内根据判决门限进行判决，实现欺骗干扰检测。本发明所公开的方法无需改变接收机结构，适用场景广泛，可应用于静态、动态接收机，对小时延欺骗干扰和复杂的拉偏式欺骗干扰具有很好的检测效果。

Description

一种基于联合SQM平方的GNSS欺骗干扰检测方法

技术领域

本发明涉及卫星信号欺骗干扰技术领域，特别涉及一种基于联合SQM(信号质量监测)平方的GNSS欺骗干扰检测方法。

背景技术

全球导航卫星系统(GNSS)已经广泛服务于各行各业，赋能各行业提质升级，深深融入到人们的日常生活和工业活动中。但GNSS设备给人们生活带来便利的同时，也带来了一定的潜在威胁。由于GNSS信号结构的开发性和信号功率的微弱性使其极易受到各种有意和无意的干扰。其中，欺骗干扰所需功率小、隐蔽性强、危害性大，受到各国军方和专家学者的广泛关注。

近些年，国内外众多高校和研究机构在欺骗干扰以及欺骗信号对接收机的影响方面做了很多研究，目前常用的欺骗干扰检测方法存在如下缺陷：(1)大多数方法改变了GNSS接收机的结构，在节省成本的条件下无法大范围推广。例如，导航信息加密的检测方法通过添加信号加密特征使得欺骗者很难获取并改变卫星的导航信息，从而增加信息安全性，但是需要国家层面对卫星信号体制进行更改，且接收机端也需要更改信号解密算法；利用欺骗干扰信号具有空间相关的特点，采用空间处理技术估计接收信号的空间特征可以实现欺骗干扰检测，但是此类方法大多需要增加接收天线数量，即便使用特殊性能的天线也需要增加整个检测系统的成本。(2)有些方法无需改变接收机的结构，但是检测适用场景较为局限，对于一些复杂精确的欺骗干扰检测效果差。例如，基于信号到达时间的检测方法对于动态性场景适用差；基于多普勒一致性的检测方法多应用于运动接收机；基于捕获阶段相关峰个数的检测方法无法应对小时延欺骗干扰；基于载噪比(C/N0)的检测方法对于逐步拉偏式的欺骗干扰检测效果差。

发明内容

为解决上述技术问题，本发明提供了一种基于联合SQM平方的GNSS欺骗干扰检测方法，该方法无需改变接收机结构，为商用接收机在节省成本的条件下实现欺骗干扰检测提供可能，且适用场景广泛，可应用于静态、动态接收机，对于小时延欺骗干扰和复杂的拉偏式欺骗干扰具有很好的检测效果。

为达到上述目的，本发明的技术方案如下：

一种基于联合SQM平方的GNSS欺骗干扰检测方法，包括如下步骤：

步骤1，获取GNSS接收机的载噪比、相关器间距、相干积分时间，并计算跟踪环路的提前路、即时路、滞后路的相干积分值；

步骤2，根据提前路、即时路、滞后路的相干积分值计算Delta、Ratio、ELP三种SQM指标；

步骤3，根据接收机的载噪比、相关器间距、相干积分时间计算SQM指标的统计特征；

步骤4，设置指标联合方式，建立联合检测量；

步骤5，设定虚警率，根据联合检测量的统计分布确定判决门限；

步骤6，设定时间窗长度，在检测时间内根据判决门限进行判决，实现欺骗干扰检测。

上述方案中，所述步骤1具体为：

A、无欺骗干扰条件下，GNSS接收机接收信号经过下变频器得到的中频信号模型为：

其中，i为卫星编号，m为卫星总数，c_i为接收到第i颗卫星的信号功率，τ_i为第i颗卫星信号的传播时延，D_i(·)为第i颗卫星调制的数据码信息，C_i(·)为第i颗卫星的C/A码序列，f_IF为中频频率，f_d,i为第i颗卫星的多普勒频移，θ_i为第i颗卫星的载波初始相位，n_fe(t)为射频前端噪声；

B、假设在相干积分时间内，接收信号的数据电平没有发生跳变，D(t-τ_i)记为采样值D(n)，n为采样时刻，且积分时间T_coh对于混频信号中的高频成分足够长，积分值近似为0，则即时路中r(t)信号经过相关器、积分滤波后的同相相干积分值可表示为：

其中，

其中，I_P,i(n)表示第i颗卫星的即时路的同相相干积分值，t₀为积分起始时刻，sinc(·)为辛格函数，T_coh为相干积分时间，

为跟踪环路对第i颗卫星的码相位跟踪值，

为跟踪环路对第i颗卫星的多普勒频率跟踪值，

为第i颗卫星的载波相位跟踪值，R(·)为C/A码的自相关函数，

为即时路的同相支路的相关噪声，

表示码相位跟踪误差；

表示频率跟踪误差；

表示相位跟踪误差；

同理，跟踪环路中其它支路的相干积分输出表示为：

其中，d为相关器间距，T_c为C/A码码片宽度，Q_P,i(n)表示第i颗卫星的即时路的正交相干积分值，I_E,i(n)表示第i颗卫星的提前路的同相相干积分值，Q_E,i(n)表示第i颗卫星的提前路的正交相干积分值，I_L,i(n)表示第i颗卫星的滞后路的同相相干积分值，Q_L,i(n)表示第i颗卫星滞后路的正交相干积分值，

为即时路的正交支路的相关噪声，

为提前路的同相支路的相关噪声，

为提前路的正交支路的相关噪声，

为滞后路的同相支路的相关噪声，

为滞后路的正交支路的相关噪声。

上述方案中，所述步骤2具体为：

根据SQM指标的定义式计算各指标数值，具体计算式如下：

其中，m_Delta,i(n)、m_Ratio,i(n)、m_ELP,i(n)分别代表第i颗卫星信号在n时刻的Delta、Ratio、ELP指标。

上述方案中，所述步骤3具体为：

A、基于步骤1的分析，针对第i颗卫星信号在n时刻的采样值，假设数据电平的采样值D(n)为1，且跟踪环路工作在稳定状态，码相位跟踪误差和载波相位误差为均0；将Δτ_i＝0，Δf_i＝0，Δθ_i＝0，D(n)＝1带入式(4)～(8)中，跟踪环路同相相干积分值和正交相干积分值的统一形式分别写为：

由于同相支路的相关噪声

与正交支路的相关噪声

相互独立，且均服从均值为0，方差为

的高斯分布，所以跟踪环路输出值的均值和方差如下所示：

E[Q_d]＝0 (15)

其中，N₀表示高斯白噪声的双边带功率谱密度，E[I_d]表示跟踪环路同相相干积分值的均值，E[Q_d]表示跟踪环路正交相干积分值的均值，D[I_d]表示表示跟踪环路同相相干积分值的方差，D[Q_d]表示跟踪环路正交相干积分值的方差，当d＜0时，I_d,Q_d分别表示提前路的同相、正交相干积分值，当d＞0时，I_d,Q_d分别表示滞后路的同相、正交相干积分值，当d＝0时，I_d,Q_d分别表示即时路的同相、正交相干积分值；

B、根据定义，SQM指标的计算描述为对不同支路的相干积分值先求比值再进行线性组合的过程，则SQM指标的一般形式表示为：

其中，

a_j表示比例系数，J表示r_j比例式的数量，x_j，y_j分别表示跟踪环路两个不同支路的相干积分输出值，假设x_j与y_j相互独立，利用泰勒公式展开和统计理论分析，计算得到r_j的统计特征为：

其中，μ_x、σ_x ²分别为变量x_j的均值和方差，μ_y、σ_y ²分别为变量y_j的均值和方差，μ_r，σ_r ²分别为变量r_j的均值和方差；

将各SQM指标变换为式(17)所示的一般形式，将式(14)～(16)带入式(19)～(20)，计算得到不同SQM指标的理论均值和方差；在高信噪比的条件下，Delta，Ratio，ELP三个SQM指标认定为服从高斯分布，均值与方差如下所示：

μ_Delta＝0 (21)

μ_ELP＝0 (25)

其中，μ_Delta表示Delta指标的均值，σ_Delta ²表示Delta指标的方差，μ_Ratio表示Ratio指标的均值，σ_Ratio ²表示Ratio指标的方差，μ_ELP表示ELP指标的均值，σ_ELP ²表示ELP指标的方差。

上述方案中，所述步骤4具体为：

利用SQM指标服从高斯分布的统计特性构建联合检测量M_cmb，其具体形式如下：

其中，M_cmb,i(n)表示第i颗卫星信号的n时刻的联合检测量，m_Delta,i(n)、m_Ratio,i(n)、m_ELP,i(n)分别代表第i颗卫星信号在n时刻的Delta、Ratio、ELP指标，μ_Delta表示Delta指标的均值，σ_Delta表示Delta指标的标准差，μ_Ratio表示Ratio指标的均值，σ_Ratio表示Ratio指标的标准差，μ_ELP表示ELP指标的均值，σ_ELP表示ELP指标的标准差；β₁表示Delta指标的系数，β₂表示Ratio指标的系数，β₃表示ELP指标的系数，β_i取值规则如下：

得到所有的M_cmb,i(n)构成的样本总体M_cmb符合以下统计特性：

M_cmb～χ²(k) (29)

其中，k＝β₁+β₂+β₃，表示参与联合检测的SQM指标数量，χ²(k)表示自由度为k的卡方分布。

上述方案中，所述步骤5具体为：

M_cmb的概率密度函数为：

式中，Γ(k/2)是Gamma函数；

M_cmb的累积分布函数为：

式中，

为不完全Gamma函数；

由此，计算虚警率为：

式中，th为判决门限，P_fa表示在无欺骗干扰攻击时，指标超过判决门限概率，即虚警率；

根据反函数存在定理，由于累积分布函数F_k(·)严格单调递增，其必存在反函数，对式(32)进行整理，得到判决门限为：

th＝F_k ^-1(1-P_fa) (33)

式中，F_k ^-1(·)为F_k(·)的反函数，设定虚警率P_fa，存在唯一的判决门限th与之对应；χ²分布对于不同的自由度和右尾概率，对应分位点的值已制成表格；通过指标联合方式确定自由度，根据需求设定虚警率，可通过查找χ²分布临界值表获得对应的判决门限。

上述方案中，所述步骤6具体为：

由于跟踪环路每隔T_coh输出一次相干积分值，所以时间窗长度设定为NT_coh，则在一次检测时间区间内存在N个联合SQM样本，P_d,i定义为第i颗卫星存在欺骗干扰攻击时联合检测量M_cmb超过门限的样本数量与样本总量的比率；

欺骗干扰检测概率P_d,i计算公式为：

其中，M_cmb,i(n)表示第i颗卫星信号的n时刻的联合检测量，th为判决门限，I(·)为指示函数，当M_cmb,i(n)＞th时，取值为1，否则，取值为0；

根据欺骗干扰检测概率P_d,i的值判断第i颗卫星是否受到欺骗干扰。

通过上述技术方案，本发明提供的一种基于联合SQM平方的GNSS欺骗干扰检测方法，有益效果如下：

1、本发明利用接收机跟踪环路中提前路，即时路，滞后路的相干积分输出值进行计算，无需改变接收机的结构，无需引入额外设备，为商用导航接收机在节省硬件成本的条件下实现欺骗干扰检测提供可能。

2、本发明解决了学者提出的基于幅度联合的联合SQM算法对欺骗干扰的检测效果相较于单一SQM指标没有改善的问题。联合SQM算法忽视了不同指标存在正偏于上门限和负偏于下门限的情况，导致联合之后异常抵消，本发明提供的方法借鉴最小二乘的思想，使用平方操作使所有的偏差变为正数再进行幅度相加，解决了上述异常抵消问题，在没有增加算法复杂度的基础上，极大提高了检测性能。

3、本发明通过识别跟踪环路相关峰的畸变判断接收机是否受到欺骗干扰攻击，对固定和运动的接收机均适用，且对于小时延欺骗干扰和复杂精确的拉偏式欺骗干扰具有很好的检测效果。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。

图1为本发明实施例所公开的一种基于联合SQM平方的GNSS欺骗干扰检测方法流程示意图。

图2是本发明实施例提供的Delta指标结果可视化图；

图3是本发明实施例提供的Ratio指标结果可视化图；

图4是本发明实施例提供的ELP指标结果可视化图；

图5是本发明实施例提供的选取Delta、Ratio、ELP三种SQM指标联合方式下的检测量M_cmb结果可视化图。

图6是本发明实施例提供的选取不同SQM指标联合方式下的欺骗干扰检测结果仿真图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。

本发明提供了一种基于联合SQM平方的GNSS欺骗干扰检测方法，如图1所示，包括如下步骤：

本实施例选取美国德州大学TEXBAT数据集中的第二条欺骗干扰数据，数据采样率为25MHz，C/A码速率为1.023MHz，欺骗信号功率恒定高于真实信号功率10dB，选取欺骗信号数据总长度为400s，欺骗攻击开始时间为第100s。使用GNSS软件接收机对数据进行处理，选取捕获结果中信号功率最高的第23号卫星，即后面公式中的i＝23，对23号卫星信号的跟踪结果进行欺骗干扰检测。

步骤1，获取GNSS接收机的载噪比、相关器间距、相干积分时间，并计算跟踪环路的提前路、即时路、滞后路的相干积分值。

A、无欺骗干扰条件下，GNSS接收机接收信号经过下变频器得到的中频信号模型为：

其中，i为卫星编号，m为卫星总数，c_i为接收到第i颗卫星的信号功率，τ_i为第i颗卫星信号的传播时延，D_i(·)为第i颗卫星调制的数据码信息，C_i(·)为第i颗卫星的C/A码序列，f_IF为中频频率，f_d,i为第i颗卫星的多普勒频移，θ_i为第i颗卫星的载波初始相位，n_fe(t)为射频前端噪声；

B、假设在相干积分时间内，接收信号的数据电平没有发生跳变，D(t-τ_i)记为采样值D(n)，n为采样时刻，且积分时间T_coh对于混频信号中的高频成分足够长，积分值近似为0，则即时路中r(t)信号经过相关器、积分滤波后的同相相干积分值可表示为：

其中，

其中，I_P,i(n)表示第i颗卫星的即时路的同相相干积分值，t₀为积分起始时刻，sinc(·)为辛格函数，T_coh为相干积分时间，

为跟踪环路对第i颗卫星的码相位跟踪值，

为跟踪环路对第i颗卫星的多普勒频率跟踪值，

为第i颗卫星的载波相位跟踪值，R(·)为C/A码的自相关函数，

为即时路的同相支路的相关噪声，

表示码相位跟踪误差；

表示频率跟踪误差；

表示相位跟踪误差；

同理，跟踪环路中其它支路的相干积分输出表示为：

其中，d为相关器间距，T_c为C/A码码片宽度，Q_P,i(n)表示第i颗卫星的即时路的正交相干积分值，I_E,i(n)表示第i颗卫星的提前路的同相相干积分值，Q_E,i(n)表示第i颗卫星的提前路的正交相干积分值，I_L,i(n)表示第i颗卫星的滞后路的同相相干积分值，Q_L,i(n)表示第i颗卫星滞后路的正交相干积分值，

为即时路的正交支路的相关噪声，

为提前路的同相支路的相关噪声，

为提前路的正交支路的相关噪声，

为滞后路的同相支路的相关噪声，

为滞后路的正交支路的相关噪声。

本实施例接收机接收信号的载噪比为53dB·HZ，相关器间距为0.5个码片，相干积分时长为1ms，相干积分值存储在维度为6×400000的跟踪矩阵中。

步骤2，根据提前路、即时路、滞后路的相干积分值计算Delta、Ratio、ELP三种SQM指标。

根据SQM指标的定义式计算各指标数值，具体计算式如下：

其中，m_Delta,i(n)、m_Ratio,i(n)、m_ELP,i(n)分别代表第i颗卫星信号在n时刻的Delta、Ratio、ELP指标。

本实施例计算得到的m_Delta,i(n)、m_Ratio,i(n)、m_ELP,i(n)可视化结果分别如图2、图3、图4所示。

步骤3，根据接收机的载噪比、相关器间距、相干积分时间计算SQM指标的统计特征。

A、基于步骤1的分析，针对第i颗卫星信号在n时刻的采样值，假设数据电平的采样值D(n)为1，且跟踪环路工作在稳定状态，码相位跟踪误差和载波相位误差为均0；将Δτ_i＝0，Δf_i＝0，Δθ_i＝0，D(n)＝1带入式(4)～(8)中，跟踪环路同相相干积分值和正交相干积分值的统一形式分别写为：

由于同相支路的相关噪声

与正交支路的相关噪声

相互独立，且均服从均值为0，方差为

的高斯分布，所以跟踪环路输出值的均值和方差如下所示：

E[Q_d]＝0 (15)

其中，N₀表示高斯白噪声的双边带功率谱密度，E[I_d]表示跟踪环路同相相干积分值的均值，E[Q_d]表示跟踪环路正交相干积分值的均值，D[I_d]表示表示跟踪环路同相相干积分值的方差，D[Q_d]表示跟踪环路正交相干积分值的方差，当d＜0时，I_d,Q_d分别表示提前路的同相、正交相干积分值，当d＞0时，I_d,Q_d分别表示滞后路的同相、正交相干积分值，当d＝0时，I_d,Q_d分别表示即时路的同相、正交相干积分值。

本实施例使用的GNSS接收机中的相关器间距为0.5个码片，当d＝-0.5时，I_d,Q_d表示提前路的同相、正交相干积分值，当d＝0.5时，I_d,Q_d表示滞后路的同相、正交相干积分值，当d＝0时，I_d,Q_d表示即时路的同相、正交相干积分值。将d的不同取值带入(14)～(16)式，计算得到相关器间距为0.5个码片时，跟踪环路中不同支路相干积分值的均值和方差，结果如下所示：

E[I_P]＝E[Q_E]＝E[Q_P]＝E[Q_L]＝0

B、根据定义，SQM指标的计算描述为对不同支路的相干积分值先求比值再进行线性组合的过程，例如，Delta指标写成

的形式，则SQM指标的一般形式表示为：

其中，

a_j表示比例系数，J表示r_j比例式的数量，x_j，y_j分别表示跟踪环路两个不同支路的相干积分输出值，假设x_j与y_j相互独立，利用泰勒公式展开和统计理论分析，计算得到r_j的统计特征为：

其中，μ_x、σ_x ²分别为变量x_j的均值和方差，μ_y、σ_y ²分别为变量y_j的均值和方差，μ_r，σ_r ²分别为变量r_j的均值和方差；

将各SQM指标变换为式(17)所示的一般形式，将式(14)～(16)带入式(19)～(20)，计算得到不同SQM指标的理论均值和方差；在高信噪比的条件下，Delta，Ratio，ELP三个SQM指标认定为服从高斯分布，均值与方差如下所示：

μ_Delta＝0 (21)

μ_Ratio＝1 (23)

μ_ELP＝0 (25)

其中，μ_Delta表示Delta指标的均值，σ_Delta ²表示Delta指标的方差，μ_Ratio表示Ratio指标的均值，σ_Ratio ²表示Ratio指标的方差，μ_ELP表示ELP指标的均值，σ_ELP ²表示ELP指标的方差。

步骤4，设置指标联合方式，建立联合检测量。

利用SQM指标服从高斯分布的统计特性构建联合检测量M_cmb，其具体形式如下：

其中，M_cmb,i(n)表示第i颗卫星信号的n时刻的联合检测量，m_Delta,i(n)、m_Ratio,i(n)、m_ELP,i(n)分别代表第i颗卫星信号在n时刻的Delta、Ratio、ELP指标，μ_Delta表示Delta指标的均值，σ_Delta表示Delta指标的标准差，μ_Ratio表示Ratio指标的均值，σ_Ratio表示Ratio指标的标准差，μ_ELP表示ELP指标的均值，σ_ELP表示ELP指标的标准差；β₁表示Delta指标的系数，β₂表示Ratio指标的系数，β₃表示ELP指标的系数，β_i取值规则如下：

得到所有的M_cmb,i(n)构成的样本总体M_cmb符合以下统计特性：

M_cmb～χ²(k) (29)

其中，k＝β₁+β₂+β₃，表示参与联合检测的SQM指标数量，χ²(k)表示自由度为k的卡方分布。β₁，β₂，β₃取值均为1时，计算得到的联合检测量M_cmb可视化结果如图5所示。

步骤5，设定虚警率，根据联合检测量的统计分布确定判决门限。

基于步骤4的分析，M_cmb服从自由度为k的χ²分布，其概率密度函数为：

式中，Γ(k/2)是Gamma函数；

M_cmb的累积分布函数为：

式中，

为不完全Gamma函数；

由此，计算虚警率为：

式中，th为判决门限，P_fa表示在无欺骗干扰攻击时，指标超过判决门限概率，即虚警率；

根据反函数存在定理，由于累积分布函数F_k(·)严格单调递增，其必存在反函数，对式(32)进行整理，得到判决门限为：

th＝F_k ^-1(1-P_fa) (33)

式中，F_k ^-1(·)为F_k(·)的反函数，设定虚警率P_fa，存在唯一的判决门限th与之对应；χ²分布对于不同的自由度和右尾概率，对应分位点的值已制成表格；通过指标联合方式确定自由度，根据需求设定虚警率，可通过查找χ²分布临界值表获得对应的判决门限。

本实施例设定虚警率为0.01，一种SQM指标参与检测，门限值为6.635；两种SQM指标参与检测，门限值为9.210；三种SQM指标参与检测，门限值为11.345。

步骤6，设定时间窗长度，在检测时间内根据判决门限进行判决，实现欺骗干扰检测。

由于跟踪环路每隔T_coh输出一次相干积分值，所以时间窗长度设定为NT_coh，则在一次检测时间区间内存在N个联合SQM样本，P_d,i定义为第i颗卫星存在欺骗干扰攻击时联合检测量M_cmb超过门限的样本数量与样本总量的比率；

欺骗干扰检测概率P_d,i计算公式为：

其中，M_cmb,i(n)表示第i颗卫星信号的n时刻的联合检测量，th为判决门限，I(·)为指示函数，当M_cmb,i(n)＞th时，取值为1，否则，取值为0；根据欺骗干扰检测概率P_d,i的值判断第i颗卫星是否受到欺骗干扰。当P_d,i大于0.5时，可认为受到欺骗干扰。

本实施例选取时间窗为100ms，在一次检测时间段内有100个联合SQM样本，采用不同的组合方式构建联合SQM检测量M_cmb进行门限判决，得到欺骗干扰检测结果如图6所示。图中横轴为接收信号时间，纵轴为欺骗干扰检测率，由于欺骗信号起初与真实信号码相位对齐，在100s到150s并未检测到欺骗存在，随时间变化，欺骗信号逐渐拉偏真实信号，150s后欺骗干扰检测概率逐渐增大，在180s时达到最大值，250s后欺骗信号与真实信号相关峰分离，欺骗干扰检测率下降为0。从图中可以看出本发明提出的联合SQM平方算法的检测概率均高于单一SQM指标的检测率，且Delta、Ratio、ELP三种SQM指标联合检测时，检测效果最好。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (7)

1.一种基于联合SQM平方的GNSS欺骗干扰检测方法，其特征在于，包括如下步骤：

步骤1，获取GNSS接收机的载噪比、相关器间距、相干积分时间，并计算跟踪环路的提前路、即时路、滞后路的相干积分值；

步骤2，根据提前路、即时路、滞后路的相干积分值计算Delta、Ratio、ELP三种SQM指标；

步骤3，根据接收机的载噪比、相关器间距、相干积分时间计算SQM指标的统计特征；

步骤4，设置指标联合方式，建立联合检测量；

步骤5，设定虚警率，根据联合检测量的统计分布确定判决门限；

步骤6，设定时间窗长度，在检测时间内根据判决门限进行判决，实现欺骗干扰检测。

2.根据权利要求1所述的一种基于联合SQM平方的GNSS欺骗干扰检测方法，其特征在于，所述步骤1具体为：

A、无欺骗干扰条件下，GNSS接收机接收信号经过下变频器得到的中频信号模型为：

其中，i为卫星编号，m为卫星总数，c_i为接收到第i颗卫星的信号功率，τ_i为第i颗卫星信号的传播时延，D_i(·)为第i颗卫星调制的数据码信息，C_i(·)为第i颗卫星的C/A码序列，f_IF为中频频率，f_d,i为第i颗卫星的多普勒频移，θ_i为第i颗卫星的载波初始相位，n_fe(t)为射频前端噪声；

B、假设在相干积分时间内，接收信号的数据电平没有发生跳变，D(t-τ_i)记为采样值D(n)，n为采样时刻，且积分时间T_coh对于混频信号中的高频成分足够长，积分值近似为0，则即时路中r(t)信号经过相关器、积分滤波后的同相相干积分值可表示为：

其中，

其中，I_P,i(n)表示第i颗卫星的即时路的同相相干积分值，t₀为积分起始时刻，sinc(·)为辛格函数，T_coh为相干积分时间，

为跟踪环路对第i颗卫星的码相位跟踪值，

为跟踪环路对第i颗卫星的多普勒频率跟踪值，

为第i颗卫星的载波相位跟踪值，R(·)为C/A码的自相关函数，

为即时路的同相支路的相关噪声，

表示码相位跟踪误差；

表示频率跟踪误差；

表示相位跟踪误差；

同理，跟踪环路中其它支路的相干积分输出表示为：

其中，d为相关器间距，T_c为C/A码码片宽度，Q_P,i(n)表示第i颗卫星的即时路的正交相干积分值，I_E,i(n)表示第i颗卫星的提前路的同相相干积分值，Q_E,i(n)表示第i颗卫星的提前路的正交相干积分值，I_L,i(n)表示第i颗卫星的滞后路的同相相干积分值，Q_L,i(n)表示第i颗卫星滞后路的正交相干积分值，

为即时路的正交支路的相关噪声，

为提前路的同相支路的相关噪声，

为提前路的正交支路的相关噪声，

为滞后路的同相支路的相关噪声，

为滞后路的正交支路的相关噪声。

3.根据权利要求2所述的一种基于联合SQM平方的GNSS欺骗干扰检测方法，其特征在于，所述步骤2具体为：

根据SQM指标的定义式计算各指标数值，具体计算式如下：

其中，m_Delta,i(n)、m_Ratio,i(n)、m_ELP,i(n)分别代表第i颗卫星信号在n时刻的Delta、Ratio、ELP指标。

4.根据权利要求3所述的一种基于联合SQM平方的GNSS欺骗干扰检测方法，其特征在于，所述步骤3具体为：

A、基于步骤1的分析，针对第i颗卫星信号在n时刻的采样值，假设数据电平的采样值D(n)为1，且跟踪环路工作在稳定状态，码相位跟踪误差和载波相位误差为均0；将Δτ_i＝0，Δf_i＝0，Δθ_i＝0，D(n)＝1带入式(4)～(8)中，跟踪环路同相相干积分值和正交相干积分值的统一形式分别写为：

由于同相支路的相关噪声

与正交支路的相关噪声

相互独立，且均服从均值为0，方差为

的高斯分布，所以跟踪环路输出值的均值和方差如下所示：

E[Q_d]＝0 (15)

其中，N₀表示高斯白噪声的双边带功率谱密度，E[I_d]表示跟踪环路同相相干积分值的均值，E[Q_d]表示跟踪环路正交相干积分值的均值，D[I_d]表示表示跟踪环路同相相干积分值的方差，D[Q_d]表示跟踪环路正交相干积分值的方差，当d＜0时，I_d,Q_d分别表示提前路的同相、正交相干积分值，当d＞0时，I_d,Q_d分别表示滞后路的同相、正交相干积分值，当d＝0时，I_d,Q_d分别表示即时路的同相、正交相干积分值；

B、根据定义，SQM指标的计算描述为对不同支路的相干积分值先求比值再进行线性组合的过程，则SQM指标的一般形式表示为：

其中，

a_j表示比例系数，J表示r_j比例式的数量，x_j，y_j分别表示跟踪环路两个不同支路的相干积分输出值，假设x_j与y_j相互独立，利用泰勒公式展开和统计理论分析，计算得到r_j的统计特征为：

其中，μ_x、σ_x ²分别为变量x_j的均值和方差，μ_y、σ_y ²分别为变量y_j的均值和方差，μ_r，σ_r ²分别为变量r_j的均值和方差；

将各SQM指标变换为式(17)所示的一般形式，将式(14)～(16)带入式(19)～(20)，计算得到不同SQM指标的理论均值和方差；在高信噪比的条件下，Delta，Ratio，ELP三个SQM指标认定为服从高斯分布，均值与方差如下所示：

μ_Delta＝0 (21)

μ_ELP＝0 (25)

其中，μ_Delta表示Delta指标的均值，σ_Delta ²表示Delta指标的方差，μ_Ratio表示Ratio指标的均值，σ_Ratio ²表示Ratio指标的方差，μ_ELP表示ELP指标的均值，σ_ELP ²表示ELP指标的方差。

5.根据权利要求1所述的一种基于联合SQM平方的GNSS欺骗干扰检测方法，其特征在于，所述步骤4具体为：

利用SQM指标服从高斯分布的统计特性构建联合检测量M_cmb，其具体形式如下：

其中，M_cmb,i(n)表示第i颗卫星信号的n时刻的联合检测量，m_Delta,i(n)、m_Ratio,i(n)、m_ELP,i(n)分别代表第i颗卫星信号在n时刻的Delta、Ratio、ELP指标，μ_Delta表示Delta指标的均值，σ_Delta表示Delta指标的标准差，μ_Ratio表示Ratio指标的均值，σ_Ratio表示Ratio指标的标准差，μ_ELP表示ELP指标的均值，σ_ELP表示ELP指标的标准差；β₁表示Delta指标的系数，β₂表示Ratio指标的系数，β₃表示ELP指标的系数，β_i取值规则如下：

得到所有的M_cmb,i(n)构成的样本总体M_cmb符合以下统计特性：

M_cmb～χ²(k) (29)

其中，k＝β₁+β₂+β₃，表示参与联合检测的SQM指标数量，χ²(k)表示自由度为k的卡方分布。

6.根据权利要求5所述的一种基于联合SQM平方的GNSS欺骗干扰检测方法，其特征在于，所述步骤5具体为：

M_cmb的概率密度函数为：

式中，Γ(k/2)是Gamma函数；

M_cmb的累积分布函数为：

式中，

为不完全Gamma函数；

由此，计算虚警率为：

式中，th为判决门限，P_fa表示在无欺骗干扰攻击时，指标超过判决门限概率，即虚警率；

根据反函数存在定理，由于累积分布函数F_k(·)严格单调递增，其必存在反函数，对式(32)进行整理，得到判决门限为：

th＝F_k ^-1(1-P_fa) (33)

式中，F_k ^-1(·)为F_k(·)的反函数，设定虚警率P_fa，存在唯一的判决门限th与之对应；χ²分布对于不同的自由度和右尾概率，对应分位点的值已制成表格；通过指标联合方式确定自由度，根据需求设定虚警率，可通过查找χ²分布临界值表获得对应的判决门限。

7.根据权利要求1所述的一种基于联合SQM平方的GNSS欺骗干扰检测方法，其特征在于，所述步骤6具体为：

由于跟踪环路每隔T_coh输出一次相干积分值，所以时间窗长度设定为NT_coh，则在一次检测时间区间内存在N个联合SQM样本，P_d,i定义为第i颗卫星存在欺骗干扰攻击时联合检测量M_cmb超过门限的样本数量与样本总量的比率；

欺骗干扰检测概率P_d,i计算公式为：

其中，M_cmb,i(n)表示第i颗卫星信号的n时刻的联合检测量，th为判决门限，I(·)为指示函数，当M_cmb,i(n)＞th时，取值为1，否则，取值为0；

根据欺骗干扰检测概率P_d,i的值判断第i颗卫星是否受到欺骗干扰。

Images