钟伦珑 刘炅坡 余亮 张卓轩

（1.中国民航大学天津市智能信号与图像处理重点实验室，天津 300300；2.中国商飞上海飞机设计研究院，上海 201210）

1 引言

考虑到全球导航卫星系统（Global Navigation Satellite System，GNSS）可以为航空提供覆盖全球的高性能导航信息，国际民航组织（ICAO）将其选作新一代区域导航系统的主用系统。但是，由于民用GNSS 信号的结构类型和扩频码字信息公开，GNSS极易受到欺骗式干扰的影响［1］，对航空安全形成威胁［2］。

按照欺骗式干扰具体表现划分，可分为突变式和缓变式两种。突变式欺骗造成的伪距突变容易被机载完好性监视感知。而缓变式欺骗在干扰初期与真实信号保持同步，当接收机跟踪欺骗信号后，采取逐步诱导的方式使伪距测量值逐渐偏离真实值，从而绕开接收机自主完好性监视。因此，对机载缓变式欺骗干扰检测技术的研究更为迫切。受机载设备适航因素等影响，机载欺骗式干扰检测技术主要分为两大类：基于信号处理的检测技术［3-5］和基于信息解算［6-8］的检测技术。基于信号处理的检测技术需改变载体接收机结构；而基于信息解算的检测技术利用检测算法，融合其他测量信息，对接收机解算输出的导航信息进行可靠性判断，可在机载区域导航计算机中实现，无需改变机载设备结构。

INS/GNSS 组合导航是区域导航的一种主用方式，可适用的基于信息解算的欺骗检测算法分为“快照式”和“连续式”两种［9］。“快照式”检测算法以当前时刻的新息作为检测量，构成检验统计量，适合检测GNSS 测量值中阶跃类型的误差，而无法检测出缓变类型的误差，例如新息欺骗检测算法和新息卡方检测算法；“连续式”检测算法将检测窗口内的新息序列作为检测量，构成检验统计量，检测缓变类型的误差，是目前常见的机载缓变式欺骗干扰检测方法，例如新息序列卡方检测算法。Wald提出的序贯概率比检测（Sequential Probability Ratio Test，SPRT）算法是一种对小偏差敏感的检测方法，广泛应用于故障检测领域中［10］，但是传统SPRT 算法检测速度较慢。为加快检测速度，Xiao［11］等提出一种分段SPRT 改进方法，减少了历史样本信息的影响，能够在短时延下正确检测故障。Wang［12］等将SPRT 算法应用到组合导航中，提出了一种联合检测算法，将SPRT 算法与模糊自适应残差检测算法相结合，在子滤波器失效的情况下，也能检测故障源并隔离故障部分，提高了系统可靠性。基于故障与欺骗产生影响的类似性，SPRT 算法可用于欺骗式干扰检测。Yuan［13］等提出一种基于SPRT 算法，将信号功率作为欺骗检测量的欺骗式干扰检测方法。但是，当欺骗信号和真实信号功率相近时，该方法检测精度变差，而且这种方法需要在信号端进行检测，改变接收机体制，适用性不好。

为了提高以新息序列为检测量的缓变式欺骗干扰检测方法的检测性能，本文受文献［13］的启发，通过将新息序列作为SPRT 算法的欺骗检测量，把SPRT 算法应用到欺骗检测当中，并提出一种自适应SPRT 算法，提高欺骗检测性能。自适应SPRT算法根据Bayes 参数估计理论，通过均衡虚警产生的风险K1和每次采样产生的风险r，利用当前时刻新息得到欺骗检验统计量的自适应补偿值，从而加快欺骗检测速度，并可通过调整风险参数K1/r进一步提高欺骗检测性能。本文提出的缓变式欺骗干扰检测算法可应用与多星受欺骗的情况，并且与现有算法相比，对伪距欺骗速率小的缓变式欺骗干扰，具有更优的检测速度和检测率。

2 缓变式GNSS欺骗干扰模型

可以通过建立基于GNSS 测量值的欺骗干扰模型，对缓变式欺骗干扰进行分析和仿真，从而避免复杂的欺骗信号仿真。载体位置对应第j颗可见卫星的伪距测量值ρGNSS，j为：

其中，cδtGNSS为GNSS 接收机钟差导致的位置偏差，δtGNSS表示GNSS 接收机时钟误差，c为光速；υρ是包括接收机内部噪声在内的GNSS 伪距测量噪声，其服从零均值的高斯分布，υρ～N（0，σG2NSS），σGNSS为伪距测量噪声标准差。根据缓变式欺骗干扰的原理，其在干扰初期通过与真实信号的码相位和载波频率保持一致，将附加的第j颗卫星的伪距欺骗量Δρj从零开始缓慢增加，缓变式GNSS 欺骗干扰模型可以表示为：

附加的第j颗卫星的伪距欺骗量Δρj可以表示为：

其中，k0为缓变式欺骗干扰在锁定接收机跟踪环路开始增加伪距欺骗量的时刻，k为当前时刻，a为伪距欺骗速率。

3 传统的新息序列卡方欺骗检测算法

卡尔曼滤波新息是实际测量向量与一步预测向量的差值［14］，可以指示测量与估计是否一致，其大小由测量误差和系统模型准确度共同决定，当系统模型准确度高时，卡尔曼滤波新息可以表征实际测量误差。历史时刻新息构成了新息序列，包含了历史测量误差［15］。新息序列卡方欺骗检测算法通过构造累加的新息序列检验统计量，相较于快照式的新息卡方检测算法，增加了对缓变式欺骗的敏感度。

设当前k时刻测量矩阵为Hk，测量向量为Zk，一步预测向量协方差矩阵为Pk，k-1，测量噪声协方差矩阵为Rk，则新息序列δZk和新息协方差矩阵Cδz，k为：

当INS/GNSS组合导航系统未受到欺骗干扰时，实际测量向量不存在偏差，由于新息序列δZk可表征实际测量误差，新息序列δZk服从高斯分布δZk～N。当INS/GNSS 组合导航系统受到欺骗干扰时，新息序列δZk服从非零均值的高斯分布δZk～N。通过检测新息序列的统计特性即可判断系统是否存在欺骗式干扰，新息序列卡方欺骗检测算法的检验统计量λk设置如下：

设可见卫星数为n，根据数理统计理论，当未受到欺骗干扰时，λk服从自由度为nk的中心卡方分布；当受到欺骗干扰时，λk服从自由度为nk的非中心卡方分布。根据卡方分布分位数的计算公式，检测门限Ts由虚警率PFA确定：

根据二元假设检验原理，判决准则为：（1）若λk＞Ts，则判断组合导航系统受到了欺骗式干扰；（2）若λk≤Ts，则判断组合导航系统正常。

4 基于自适应SPRT的欺骗检测算法

4.1 新息序列SPRT欺骗检测算法

新息序列卡方欺骗检测算法只能判断组合导航系统是否受到欺骗，但无法识别具体哪路卫星信号受到欺骗，并且对小伪距欺骗速率的缓变式欺骗的检测性能不高［16］。SPRT 算法可依次对可见卫星进行序贯似然比检测，基于新息序列，本文提出一种基于SPRT 算法的欺骗检测算法，将新息序列δZk的第j行作为检测量对第j颗卫星是否受到欺骗做出检测，具体方法如下：

设δZk，（jj=1，2，…，n）为新息序列δZk的第j行，表示当前k时刻的第j颗卫星的新息序列。历史时刻到当前k时刻的新息序列δZk，j构成一组样本{ft，j=δZt，j|t=1，2，…，k｝。卫星受欺骗会造成伪距偏差的存在，导致伪距差量测值的均值改变，但其依然服从高斯分布，这组样本的样本均值fˉk，j和方差的当前值为：

基于二元假设检验［8］，定义原假设H0：第j颗可见卫星未受欺骗，ft，j(t=k)=0；备择假设H1：j颗可见卫星受欺骗，ft，j(t=k)=fˉk，j。在二元假设的条件下，观测样本ft，j的条件概率密度函数可表示为：

欺骗检验统计量的当前值λk，j可表示为历史时刻到当前k时刻的对数似然比Zt，（jt=1，2，…，k）之和，即：

新息序列SPRT 欺骗检测算法的检测门限Tρ由虚警率α和漏警率β计算得到［17］：

将计算得到的第j颗可见卫星欺骗检验统计量λk，j与欺骗检测门限Tρ进行比较，并根据比较结果判断第j颗可见卫星是否受到欺骗，判决准则为：（1）若λk，j＞Tρ，则判断第j颗可见卫星受到欺骗，并向系统告警；（2）若λk，j≤Tρ，则判断第j颗可见卫星正常，继续采样。

4.2 自适应SPRT算法

新息序列SPRT 算法可以对每路卫星是否受到欺骗做出检测。但是，当伪距欺骗量较小时，由式（11）可知，SPRT 欺骗检验统计量λk，j需较长时间才能大于欺骗检测门限，导致检测时间过长，难以满足飞行阶段机载导航的告警时间需求。针对这一问题，本文应用Bayes 参数估计理论，提出一种自适应SPRT 改进算法。综合考虑虚警风险、每次采样所付出的风险和当前k时刻欺骗检测量的大小，在Bayes 平均风险最小条件下，计算式（11）的SPRT欺骗检验统计量的自适应补偿值，加快欺骗检测速度。

设原假设H0的先验概率为π0；备择假设H1的先验概率为π（1π0+π1=1）；每次采样所付出的风险为r；H1为真但接受H0所产生的漏警风险为K0；H0为真但接受H1所产生的虚警风险为K1；虚警率为α；漏警率为β；H0为真，预期的采样次数为N0；H1为真，预期的采样次数为N1；根据Bayes参数估计理论［18］，当H0为真时，所付出的Bayes平均风险R0为：

当H1为真时，所付出的Bayes平均风险R1为

由上两式，总Bayes平均风险R为：

在传统SPRT 算法中，预期采样次数N0和由预设的虚警率α和漏警率β、当前k时刻的对数似然比的期望值EH0(Zk，j)和EH1(Zk，j)计算得来：

此外，可以推导出传统SPRT 算法的上门限lnA、下门限lnB［20］和虚警率α、漏警率β的关系式为：

按照机载GNSS 接收机性能要求规范［21］，虚警率α和漏警率β规定值接近于0，且α＜＜β，故式（17）、（18）、（19）可近似为：

将式（20）代入式（16）可得当前k时刻总的Bayes平均风险R为：

式（21）相对于A求偏导，并令其偏导为零，即可得到使总Bayes平均风险最小的Ak：

可以认为第j颗可见卫星是否受欺骗干扰的先验概率是相等的，即π0=π1，代入式（22），将使总Bayes 平均风险最小的Ak的对数lnAk作为自适应补偿值ζk，j：

其中，风险参数K1/r设置较小时，在相同的欺骗量下会做出更小的补偿；K1/r设置较大时，在相同的欺骗量下会做出更大的补偿。风险参数K1/r代表对欺骗干扰施加欺骗量的敏感程度，设置K1/r=0.2。EH1(Zk，j)为第j颗可见卫星受欺骗时对数似然比的期望值：

EH1(Zk，j)的大小随欺骗干扰加入的伪距欺骗量变化，因此补偿值ζk，j的大小也将随伪距欺骗量而变化，可实现对检验统计量λk，j的自适应补偿，减小达到检测门限的时间。在各个采样时刻，由式（23）补偿式（11）的统计量，可得基于本文改进SPRT 算法的第j颗可见卫星欺骗检验统计量

写成递推形式为：

当多路卫星信号受到低速率的缓变式欺骗时，由于缓变式欺骗影响下的GNSS 测量值误差随时间累积的快慢与伪距欺骗速率有关，伪距欺骗速率越小，误差累积的越慢，新息变化的就越慢，欺骗检验统计量到达门限值的时间就越久，最终导致欺骗检测性能降低。因此，多路卫星信号受到低速率的缓变式欺骗时，欺骗检测性能会显著降低。此时，可通过调节风险参数K1/r，增加对伪距欺骗量的敏感度，在同样的伪距欺骗量下做出更大的补偿，提高在多星受欺骗时的检测性能。

本文提出的基于自适应SPRT 的欺骗检测算法的算法流程可总结如图1所示。

图1 基于自适应SPRT的欺骗检测算法流程图Fig.1 Flowchart of spoofing detection algorithm based on adaptive SPRT

首先根据4.1 所提新息序列SPRT 欺骗检测算法计算SPRT欺骗检验统计量；然后根据4.2所提自适应SPRT 算法计算自适应补偿值，补偿SPRT 欺骗检验统计量构成自适应SPRT 欺骗检验统计量；接着与欺骗检测门限Tρ对比，判断第j颗卫星是否受到欺骗；最后遍历所有可见卫星，完成对当前k时刻所有可见卫星的检测。

5 仿真实验

针对INS/GNSS紧耦合模式进行仿真实验，设计以下了3个实验：

1）在单颗卫星受到伪距欺骗速率为0.1 m/s、0.2 m/s、0.3 m/s、0.4 m/s 的缓变式欺骗时比较新息序列卡方欺骗检测算法、新息序列SPRT 欺骗检测算法和自适应SPRT 欺骗检测算法的检测时间，验证本文自适应SPRT 算法在检测速度上的优势。并且以单星受到伪距欺骗速率为0.1 m/s 的缓变式欺骗为例，利用自适应SPRT 欺骗检测算法识别受欺骗卫星。

2）在单颗卫星受到缓变式欺骗时比较上述三种算法的欺骗检测率，验证本文自适应SPRT 算法在检测率上的优势。

3）在多颗卫星受到缓变式欺骗时，改变参数K1/r的取值，验证参数K1/r对本文自适应SPRT 算法检测性能的影响。

5.1 仿真设置

根据前面理论分析，将新息序列卡方欺骗检测算法、新息序列SPRT 欺骗检测算法以及本文自适应SPRT 欺骗检测算法嵌入到INS/GNSS 组合导航系统仿真平台中。模拟巡航阶段受欺骗干扰，进行蒙特卡洛仿真，分析欺骗式干扰检测性能。

设置30 颗GNSS 卫星均匀分布在6 个等圆轨道上。根据ICAO 附件10［22］，GNSS 模块参数设置如表1 所示；IMU 模块参数设置如表2 所示；卡尔曼滤波器参数设置如表3所示。

表1 GNSS模块参数Tab.1 Parameters of GNSS model

表2 IMU模块参数Tab.2 Parameters of IMU model

表3 卡尔曼滤波器参数Tab.3 Parameters of Kalman filter

设置遮蔽角为10°，在模拟飞机巡航阶段航迹中，保持可见卫星数为8颗，并设置不同的伪距欺骗速率。欺骗场景设置如表4所示。

表4 欺骗场景设置Tab.4 Spoofing scene setup

5.2 评价指标

为了综合反映欺骗检测速度和检测精度，表征算法综合检测性能，设置计算欺骗检测率Pd，其计算公式为：

其中，q为蒙特卡洛实验次数，td为检测时间，表示成功检测的次数。参照ICAO 附件10 中对机载设备巡航阶段30 秒告警时间的要求，将加入欺骗式干扰的30秒内检验统计量超过欺骗检测门限的记为一次成功检测。

5.3 仿真结果与分析

1）实验1

对在200 s 开始对卫星PRN1 施加伪距欺骗速率为0.1 m/s、0.2 m/s、0.3 m/s、0.4 m/s 的缓变式欺骗进行仿真实验，验证本文算法在检测速度上的优越性。三种算法的欺骗检验统计量随时间的变化曲线图如图2、图3、图4和图5所示。

图2、图3、图4和图5中的纵坐标代表三种检测算法的检验统计量，图中新息序列卡方检测算法、新息序列SPRT 算法和自适应SPRT 算法的检验统计量分别由公式（6）、（11）和（25）所计算。通过图2、图3、图4 和图5 可以看出，随着伪距欺骗速率增加，三种方法的检验统计量的增加速度都变快了，超过检测门限所需的时间减小。新息序列卡方检测算法随着时间的增加，其检验统计量所服从的卡方分布的自由度增加，导致其检测门限会随时间增大。这与其随时间累加的检验统计量相符合，当欺骗存在时，其检验统计量的累加速度会超过检测门限随时间增大的速度，在欺骗存在一段时间后其检验统计量会超过检测门限。这种累加的方法可以累计微小的欺骗量，实现对小欺骗量的检测，但这也导致其超过检测门限的时间长。自适应SPRT算法与新息序列SPRT 算法相比，自适应SPRT 算法在新息序列SPRT 算法的基础上根据欺骗量的大小做出自适应补偿，在相同的欺骗环境下，自适应SPRT算法的检验统计量能更快地超过检测门限。

图2 伪距欺骗速率0.1 m/sFig.2 Pseudorange spoofing rate 0.1 m/s

图3 伪距欺骗速率0.2 m/sFig.3 Pseudorange spoofing rate 0.2 m/s

在第200 s 加入欺骗之前，本文改进SPRT 欺骗检验统计量存在小幅波动，这是由于伪距测量噪声的存在造成的，自适应补偿值因噪声功率小，不会形成较大补偿产生虚警。在200 s 加入欺骗后，缓变式欺骗初始加入的欺骗量小，欺骗检验统计量与未加入欺骗时相比不会有太大变化，随着欺骗量的增加欺骗检验统计量会逐渐超过检测门限，但由于三种算法的欺骗检验统计量计算方式不同造成了欺骗检验统计量超过检测门限的所需的时间不同。

通过图2、图3、图4 和图5，可列出0.1 m/s、0.2 m/s、0.3 m/s、0.4 m/s 不同伪距欺骗速率时，三种算法的欺骗检测时间，如表5所示。

图4 伪距欺骗速率0.3 m/sFig.4 Pseudorange spoofing rate 0.3 m/s

图5 伪距欺骗速率0.4 m/sFig.5 Pseudorange spoofing rate 0.4 m/s

通过表5 可以看出，施加0.1 m/s、0.2 m/s、0.3 m/s、0.4 m/s 不同伪距欺骗速率时，在欺骗检测速度方面，本文自适应SPRT 算法优于新息序列SPRT 算法，优于新息序列卡方检测算法。本文自适应SPRT 算法对小欺骗量更加敏感，在检测速度上具有优越性。

表5 三种算法的欺骗检测时间Tab.5 Spoofing detection time of the three algorithms

在单星受到伪距欺骗速率为0.1 m/s 的缓变式欺骗情况下，全部8颗可见卫星的自适应SPRT欺骗检验统计量随时间变化曲线图如图6所示。

通过图6 可以看出，在200 s 加入欺骗后，可见卫星PRN1 的检验统计量超过检测门限，而其他可见卫星的检验统计量会有小幅度的增加，但并未超过检测门限，可以识别出受欺骗的卫星为PRN1。在可见卫星PRN1 受到欺骗，但其他可见卫星的检验统计量会小幅度增加的原因在于，单颗星受到的欺骗所引入的欺骗量会通过滤波器解算过程，使状态变量缓慢偏离真实值，从而致使未受到欺骗卫星的新息序列检测量偏离真实值，最终导致未受到欺骗卫星的检验统计量小幅度增加，但这种变化与受欺骗卫星相比较微小，不影响检测算法对受欺骗卫星的识别。

图6 单星受欺骗时自适应SPRT欺骗检验统计量随时间变化曲线Fig.6 Variation curve of adaptive SPRT spoofing test statistics when single star is spoofed

2）实验2

在200 s 开始对卫星PRN1 施加不同伪距欺骗速率的缓变式欺骗，进行2000 次蒙特卡洛实验，由式（27）计算欺骗检测率。通过比较30 s 告警时间约束下的欺骗检测率，验证本文算法在检测精度和检测速度上的综合性能优势。新息序列卡方欺骗检测算法、新息序列SPRT 欺骗检测算法和本文自适应SPRT 欺骗检测算法在单星受欺骗时欺骗检测率随伪距欺骗速率的变化曲线如图7所示。

图7 三种算法的欺骗检测率变化曲线Fig.7 Variation curve of spoofing detection rate of the three algorithms

通过图7可以看出，本文自适应SPRT算法在较小的0.1 m/s 到0.7m/s 伪距欺骗速率下的欺骗检测率明显优于新息序列SPRT 算法和新息序列卡方检测算法，0.4 m/s伪距欺骗速率时本文方法的欺骗检测率已经收敛到100%，而新息序列SPRT 算法和新息序列卡方检测算法分别在0.9 m/s和1.4 m/s伪距欺骗速率时欺骗检测率收敛到100%。对于较小伪距欺骗速率的缓变式欺骗，与新息序列SPRT 算法和新息序列卡方检测算法相比，本文方法依然能够有效地检测，并且欺骗检测率的收敛速度也要快于其他两种算法。验证了本文自适应SPRT 算法在检测精度和检测速度上的综合性能优势。

3）实验3

通过给多颗卫星施加缓变式欺骗，并调节K1/r，进行2000次蒙特卡洛实验计算欺骗检测率，验证在多颗卫星受到缓变式欺骗时，参数K1/r的取值对自适应SPRT 欺骗检测算法检测性能的影响。本文自适应SPRT 欺骗检测算法在多星受欺骗时不同风险参数下，欺骗检测率随伪距欺骗速率的变化曲线如图8所示。

图8 多星受欺骗时不同参数下欺骗检测率的变化曲线Fig.8 Variation curve of spoofing detection rate under different parameters when multiple satellites are spoofed

通过图7和图8可以看出，在风险参数K1/r设置为原始值0.2，相比于单星受欺骗，多星受欺骗情况下的欺骗检测率更低。单星受欺骗时本文算法在0.4 m/s伪距欺骗速率时，欺骗检测率收敛到100%；而在多星受欺骗时，在0.8 m/s 伪距欺骗速率时，欺骗检测率才收敛到100%。

通过图8可以看出，在多星受欺骗场景下，设置更大的参数K1/r=0.5，K1/r=1.0，K1/r=2.0，本文自适应SPRT 算法的欺骗检测率有所增加。通过调整参数K1/r，算法对欺骗量的敏感程度更高，提高了在多星受欺骗时的欺骗检测率。但是由于各类噪声的存在，若持续增大参数K1/r，过于敏感的自适应补偿值会做出过大的补偿，检测算法可能会将噪声误判为欺骗式干扰，导致虚警，所以参数K1/r的选定需要视应用中实际情况而定。

6 结论

本文利用新息序列，结合Bayes 参数估计理论，提出了一种适用于多星受欺骗场景并对小伪距欺骗量敏感的缓变式欺骗干扰检测方法。实验结果表明，本文方法在保证检测速度满足告警时间要求的前提下，在多星受欺骗和伪距欺骗速率较小情况均具有较高的欺骗检测率。与常见的INS/GNSS 组合导航系统欺骗检测算法相比，对于隐蔽性强的缓变式欺骗，本文方法体现出了更好的检测性能。