GNSS欺骗攻击指利用专门的欺骗设备播发与真实信号类似的欺骗信号，使目标接收机不再跟踪真实卫星信号，输出错误的伪距观测量，进而诱导目标接收机解算结果偏离真实位置^[1]。相较其他无线电干扰，欺骗干扰具有隐蔽性高，危害大的特点。因此，改善接收机的抗欺骗能力，是提高卫星导航可靠性与安全性的重要技术。

目前，常见的抗欺骗技术主要分为3类。第1类抗欺骗技术主要基于天线技术，包括双天线抗欺骗技术^[2-4]、多天线抗欺骗技术^[5-9]、以及移动单天线抗欺骗技术^[10]。该类技术针对不同欺骗场景的检测成功率均保持于较高水平且场景适应能力强，然而劣势在于令接收机增加天线数目或相应伺服控制机构，提高硬件成本，普遍适用性较差。第2类抗欺骗技术主要基于观测量的一致性，如基于相关峰畸变的欺骗监测技术^[11-12]，基于伪距观测量一致性的RAIM(Receiver Autonomous Integrity Monitoring)的抗欺骗技术^[13-14]、基于多普勒频率与伪距变化率一致性的欺骗检测技术^[15]，基于钟差变化的欺骗检测技术^[16]、基于多传感器敏感导航信息一致性的组合导航抗欺骗技术^[17-19]等。该类技术通过伪距、钟差、多普勒等观测量之间的一致性差异检测GNSS定位解算异常，但往往局限于欺骗监测而难以获取真实的定位解算结果。第3类抗欺骗方法主要基于基带改进方案，依据欺骗信号与多径信号具有一定的相似性，将一些抗多径算法试用在抗欺骗领域，包括MEDLL(Multipath Estimation Delay Locked Loop)^[20-21]，CADLL(Coupled Amplitude Delay Lock Loop)^[22-23]、SML(Sequential Maximum Likelihood)^[24]及其改进MSML(Modified Sequential Maximum Likelihood)^[25]，同时估计真实信号与欺骗信号码延时与幅值，再利用某些假设条件，如欺骗信号码延时大于真实信号，剔除欺骗信号参数，根据余下的真实码延时获得真实伪距，最后估计出真实位置。该类方法的核心在于信号参数的估计与真实信号的辨识。

参数估计方法往往需要大量相关器拟合实际相关曲线，为了减少相关器的数量，本文提出了一种基于卡尔曼滤波器的跟踪环路，利用九路相关器即可同时估计真实信号与欺骗信号的相关幅值与码延时。现有参数估计方法在抗多径影响时往往假设直射信号强度最大且首先到达接收机，而在欺骗干扰下该假设条件不一定成立，对此本文提出基于改进RAIM的辨识方法，最终获得一组真实伪距解算接收机真实位置，以抑制欺骗干扰对接收机定位结果的诱导效果，增强接收机在欺骗干扰下的可用性。在所有卫星信号均被欺骗的极端条件下，该方法由于可以保持真实信号参数的估计和真实伪距的输出，使得最终定位结果仍可回到真实位置，避免接收机定位结果被欺骗干扰影响。

1 欺骗干扰对GNSS接收机相干积分值影响分析

欺骗信号与GNSS信号具有类似的信号结构和伪码序列，但具有不同的码延时和信号功率。因此，欺骗干扰下，不考虑欺骗信号对载波多普勒频率的影响与噪声项，以GPS(Global Positioning System)载波L1民用C/A码为例，接收机接收到一颗卫星的中频信号可写为^[26]

$ \begin{array}{*{20}{l}} {s(t) = {s^{{\rm{au}}}}(t) + {s^{{\rm{sp}}}}(t) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sqrt {2{P^{{\rm{au}}}}} C(t - {\tau ^{{\rm{au}}}}(t))D(t - {\tau ^{{\rm{au}}}}(t)) \cdot }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{sin}}[2\pi {f^{{\rm{au}}}}(t - {\tau ^{{\rm{au}}}}(t))] + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sqrt {2{P^{{\rm{sp}}}}} C(t - {\tau ^{{\rm{sp}}}}(t))D(t - {\tau ^{{\rm{sp}}}}(t)) \cdot }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{sin}}[2\pi {f^{{\rm{sp}}}}(t - {\tau ^{{\rm{sp}}}}(t))]} \end{array} $

（1）

式中：上标au代表GNSS真实信号相关参数；上标sp代表欺骗信号相关参数；P为信号功率；τ为信号码延时；C为C/A码序列；D为导航数据；f为接收信号的载波频率经接收机前端下变频后的频率，由标称中频频率与多普勒频率组成。

在欺骗信号的影响下，接收机如能保持跟踪，则产生一组正交的本地信号$u_{{\rm{os}}}^*\left( t \right) = \sqrt 2 C\left( {{\delta _i} + t - \tau \left( t \right)} \right){\rm{exp}}\left[ {{\rm{j}}2{\rm{ \mathsf{ π} }}{f^0}\left( {t - \tau \left( t \right)} \right)} \right] $与接收信号混合，即u_f(t)=s(t)u_os^*(t)，其中，δ_i表示第i路本地码延时，传统接收机一般采用即时、滞后、提前3路码延时，取值分别为0、0.5、-0.5码片；τ(t)表示本地复制码延时；f⁰为本地复制载波频率。并经过积分，得到同相I支路与正交Q支路相干积分值。在欺骗干扰下，第i路相干积分值可写为

$ \begin{array}{*{20}{l}} {{I^{{\delta _i}}} = {\rm{Re}} (\int_0^{{T_{{\rm{ coh }}}}} {{u_{\rm{f}}}} (t){\rm{d}}t) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A^{{\rm{ au }}}}R({\delta _i} + {\alpha ^{{\rm{au}}}} - \Delta {\tau ^{\rm{s}}}){\rm{cos}}(\Delta \varphi ) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A^{{\rm{sp}}}}R({\delta _i} + {\alpha ^{{\rm{sp}}}} - \Delta {\tau ^{\rm{s}}}){\rm{cos}}(\Delta \varphi )} \end{array} $

（2）

$ \begin{array}{*{20}{l}} {{Q^{{\delta _i}}} = {\rm{Im}} (\int_0^{{T_{{\rm{ coh }}}}} {{u_{\rm{f}}}} (t){\rm{d}}t) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A^{{\rm{au}}}}R({\delta _i} + {\alpha ^{{\rm{au}}}} - \Delta {\tau ^s}){\rm{sin}}(\Delta \varphi ) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {A^{{\rm{sp}}}}R({\delta _i} + {\alpha ^{{\rm{sp}}}} - \Delta {\tau ^s}){\rm{sin}}(\Delta \varphi )} \end{array} $

（3）

式中：I^δ_i与Q^δ_i为第i路同相与正交相相干积分值；Re(·)与Im(·)分别为取实部与取虚部；T_coh为相干积分时间；$ A = \sqrt P {\rm{sinc}}({\omega _{\rm{e}}}{T_{{\rm{coh}}}}/2)$为自相关最大幅值，ω_e为本地复制载波与接收机信号载波角频率差；R(·)为伪码自相关函数；α^au为在不考虑跟踪误差时接收机本地复现即时码码延时与真实信号码延时之差，即α^au=τ^au-τ；α^sp为在不考虑跟踪误差时接收机本地复现即时码码延时与欺骗信号码延时之差，即α^sp=τ^sp-τ；Δτ^s为跟踪误差，表示接收机实际接收信号码延时与本地即时码的相位差；Δφ为载波相位误差，当欺骗信号与真实信号载波完全一致，且接收机稳定跟踪时其值近似为0，否则该值增加，对应同相相干积分值I^δ_i减小而正交相相干积分值Q^δ_i增大，在估计信号幅值时容易被自相关最大幅值项吸收，影响幅值估计结果，但本身与码延时参数不相关，影响相对较少。由于积分时间内避免导航数据反相，即导航数据在一次相干积分时间内不会跳变，仅影响积分结果的符号，为使公式更简洁，此处省略导航数据项。

式(2)和式(3)所示的同相与正交相相干积分值通过码鉴相器可获得码跟踪误差，常用码鉴相器与E-L相关，其中$ E = \sqrt {{{({I^{ - 0.5}})}^2} + {{({Q^{ - 0.5}})}^2}} , L = \sqrt {{{({I^{0.5}})}^2} + {{({Q^{0.5}})}^2}} $，上标代表超前、滞后支路相关器的相对码延时。鉴相器测量的码跟踪误差再经过一阶或二阶低通滤波器，滤除噪声后控制伪码生成器生成本地码。整个跟踪流程形成反馈环路，待跟踪环路稳定后，鉴相器输出跟踪误差Δτ^s∝E-L→0。

在无欺骗干扰时，接收信号仅为真实信号，跟踪环路稳定后，即完全跟踪上真实信号，此时跟踪误差Δτ^s=τ^au-τ≈0；在欺骗攻击突然出现的瞬间，接收机由于前一时刻仍跟踪在真实信号上，此时α^au=0，但由于欺骗干扰的突然出现，相关值瞬间变化，E与L不再近似相等，Δτ^s会有瞬时的增大。但该状态不会持续太久，反馈跟踪环路会由于误差的增加而不断自我调节，使Δτ^s逐渐减小，重新达到新的跟踪稳定状态，因此，接收机实际跟踪点会因欺骗信号影响而发生偏移^[27]，当跟踪误差Δτ^s再次趋于0时，接收机实际跟踪的是真实信号与欺骗信号的混合效果，该混合效果不仅与欺骗信号与真实信号的相对延时相关，还与两信号的幅值比相关。

图 1展示了欺骗干扰下接收机跟踪环路重新进入稳态后，实际跟踪的码延时与真实信号、欺骗信号码延时之间的关系。以传统的3路本地码为例，忽略跟踪误差Δτ^s(即Δτ^s=0)，接收机重新进入稳态，码延时0.5 chip处相干积分值L与-0.5 chip处相干积分值E必然相等，即E-L=0，此时，接收机跟踪点在码延时0 chip所在处，未与真实信号峰对齐，也并未与欺骗峰对齐。

真实信号自相关峰(真实位置)与欺骗信号自相关之间的码延时间隔α为欺骗系统预设的相对欺骗码延时距离，有α=α^au+α^sp。当该值改变时，混合信号对应的相关函数也会发生变化，因此接收机具体跟踪点会受α影响。类似地，真实信号与欺骗信号的自相关幅值A^au与A^sp也会影响混合相关，从而对最终的跟踪点造成影响。因此，在欺骗干扰下，相关幅值与码延时皆是必要参数，需要同时进行估计。

2 基于EKF的真实与欺骗信号参数估计环路

从图 1的示意中可以看出，相关函数的形状对欺骗信号较敏感，而鉴相器在欺骗干扰下会迅速失效，不能长时间感知相关函数畸变，因此可以直接利用多相关器输出估计真实信号与欺骗信号参数。然而相关器运算量巨大，因此为了尽可能减少相关器数量，利用EKF(Extended Kalman Filter)进行信号参数估计，估计结果经低通滤波器降低噪声影响，输入伪码与载波NCO(Number Controlled Oscillator)，参数估计环路结构如图 2所示。

由于欺骗干扰下接收机接收到真实信号与欺骗信号的混合信号相干积分结果不仅与码延时误差有关，还与真实信号、欺骗信号各自的幅值和本地跟踪点之间的码延时有关。因此，基于EKF的欺骗信号参数估计算法包含6个状态量，2维用于码跟踪的码参数，即本地即时码与实际欺骗位置的码延时误差Δτ^s与变化率$ \Delta {{\dot \tau }^{\rm{s}}}$；真实信号与欺骗信号幅值与相对延时参数各2维，即真实信号幅值A^au、欺骗信号幅值A^sp、实际欺骗相对码延时α^au与实际欺骗偏差α^sp，因此状态方程的状态量X表示为

$ \mathit{\boldsymbol{X}} = {\left[ {\begin{array}{*{20}{l}} {\Delta {\tau ^{\rm{s}}}}&{\Delta {t^{\rm{s}}}}&{{A^{{\rm{au}}}}}&{{A^{{\rm{sp}}}}}&{{\alpha ^{{\rm{au}}}}}&{{\alpha ^{{\rm{sp}}}}} \end{array}} \right]^{\rm{T}}} $	（4）

进一步状态方程可写作

$ {\mathit{\boldsymbol{X}}_{k + 1}} = {\mathit{\boldsymbol{F}}_k}{\mathit{\boldsymbol{X}}_k} + {\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{U}}_k} + {\mathit{\boldsymbol{W}}_k} $	（5）

式中：F_k为系统状态量在k历元到k+1历元的一步状态转移矩阵，为X提供状态更新的条件；W_k为9×1阶的过程噪声向量；U_k为输入控制量；B_k为与状态量对应的输入关系矩阵。其中状态转移矩阵F_k和输入关系矩阵B_k写作

$ {\mathit{\boldsymbol{F}}_k} = \left[ {\begin{array}{*{20}{l}} 1&T&0&0&0&0\\ 0&1&0&0&0&0\\ 0&0&1&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&1&0&0\\ 0&1&0&0&0&1 \end{array}} \right] $	（6）

$ {\mathit{\boldsymbol{B}}_k} = {\left[ {\begin{array}{*{20}{l}} { - T}&0&0&0&0&0 \end{array}} \right]^{\rm{T}}} $	（7）

考虑本地即时码与实际欺骗位置的码延时误差Δτ_k^s对相关函数的实际影响，建立9维观测方程：

$ {\mathit{\boldsymbol{Z}}_k} = \left[ {\begin{array}{*{20}{l}} {[A_k^{{\rm{au}}}R({\delta _1} + a_k^{{\rm{au}}} - \Delta \tau _k^{\rm{s}}) + A_k^{{\rm{sp}}}R({\delta _1} + a_k^{{\rm{sp}}} - \Delta \tau _k^{\rm{s}})]{\rm{cos}}(\Delta {\varphi _k})}\\ {[A_k^{{\rm{au}}}R({\delta _1} + a_k^{{\rm{au}}} - \Delta \tau _k^{\rm{s}}) + A_k^{{\rm{sp}}}R({\delta _1} + a_k^{{\rm{sp}}} - \Delta \tau _k^{\rm{s}})]{\rm{sin}}(\Delta {\varphi _k})}\\ {[A_k^{{\rm{au}}}R({\delta _2} + a_k^{{\rm{au}}} - \Delta \tau _k^{\rm{s}}) + A_k^{{\rm{sp}}}R({\delta _2} + a_k^{{\rm{sp}}} - \Delta \tau _k^{\rm{s}})]{\rm{cos}}(\Delta {\varphi _k})}\\ {[A_k^{{\rm{au}}}R({\delta _2} + a_k^{{\rm{au}}} - \Delta \tau _k^{\rm{s}}) + A_k^{{\rm{sp}}}R({\delta _2} + a_k^{{\rm{sp}}} - \Delta \tau _k^{\rm{s}})]{\rm{sin}}(\Delta {\varphi _k})}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \vdots }\\ {[A_k^{{\rm{au}}}R({\delta _i} + a_k^{{\rm{au}}} - \Delta \tau _k^{\rm{s}}) + A_k^{{\rm{sp}}}R({\delta _i} + a_k^{{\rm{sp}}} - \Delta \tau _k^{\rm{s}})]{\rm{cos}}(\Delta {\varphi _k})}\\ {[A_k^{{\rm{au}}}R({\delta _i} + a_k^{{\rm{au}}} - \Delta \tau _k^{\rm{s}}) + A_k^{{\rm{sp}}}R({\delta _i} + a_k^{{\rm{sp}}} - \Delta \tau _k^{\rm{s}})]{\rm{sin}}(\Delta {\varphi _k})} \end{array}} \right] + \mathit{\boldsymbol{V}}_k^{{\rm{es}}} $

（8）

式中：V_k^es表示9×1阶的观测噪声。第j=2, 5, 8路相关器为主相关器，码延时δ₂、δ₅、δ₈分别为-0.5、0、0.5 chip。其他相关器为辅助相关器，为提高本地码序列生成效率，考虑接收机前端采样率f_s=16.367667 MHz，相关器间距选择1.023 Mbps/f_s≈0.0625 chip的整数倍，因此δ₁、δ₄、δ₇较主相关器提前0.0625 chip，对应码延时为-0.5625、-0.0625、0.4375 chip，δ₃、δ₆、δ₉则较主相关器滞后0.0625 chip，为-0.4375、0.0625、0.5625 chip，主相关器对应观测点与各信号成分之间的相位关系如图 3所示。

该跟踪环路实际跟踪信号是真实信号与欺骗信号的混合信号，直接输出的伪距既非真实伪距，也非欺骗伪距，因此利用状态方程估计的α^au与α^sp，可以同时估计真实信号与欺骗信号的码延时，进而获得真实伪距与欺骗伪距，即

$ \left\{ {\begin{array}{*{20}{l}} {\rho _i^1 = {\rho _i} - \gamma _\alpha ^{{\rm{au}}}}\\ {\rho _i^2 = {\rho _i} + \gamma _\alpha ^{{\rm{sp}}}} \end{array}} \right. $	（9）

式中：ρ_i为受欺骗影响通道解算伪距；γ为码延时与伪距之间的关系系数，对于C/A码有γ=293.3 m/chip。所获得的两伪距差正比于相对码延时，即ρ_i²-ρ_i¹=γα。需要指出的是，跟踪环路本身估计的α^au与α^sp并非正确对应真实信号与欺骗信号相对跟踪点的码间隔。因此，式(9)只能获得2组伪距，并不能直接分辨出哪一组是真实伪距，因此需要进一步对获得的两组伪距进行辨识。

3 基于最小残差的欺骗伪距辨识算法

辨识算法与基于RAIM的欺骗检测算法类似，但由于受欺骗卫星的跟踪通道将输出2组伪距，一组为真实伪距，另一组必为欺骗伪距，因此针对该特点对经典基于RAIM的欺骗检测算法进行改进。

基于最小二乘残差的RAIM算法以包含残余平方加权和的标准检测量来监测各个通道测量值之间的一致性程度^[28]。而一致性检测要求观测量数量多于解的维数，即可见星数目不低于5颗。对于N个受欺骗影响通道与M个未受欺骗影响通道(N+M≥5)：受欺骗干扰通道输出两个伪距，其对应跟踪通道输出两个伪距；不受欺骗干扰通道输出一个伪距。伪距组合观测量形式为

$ {\rho _i} = \left\{ {\begin{array}{*{20}{l}} {\rho _i^1,\rho _i^2}&{i = 1,2, \cdots ,N}\\ {{\rho _i}}&{i = N + 1,N + 2, \cdots ,N + M} \end{array}} \right. $	（10）

最终可获得2^N种组合，具体如表 1所示。

对于第l组伪距组合，可获得一组伪距残差w_b^com_l，进一步计算辨识统计量：

$ T_{{\rm{dis}}}^{{\rm{co}}{{\rm{m}}_l}} = \sqrt {{{(\mathit{\boldsymbol{w}}_{\rm{b}}^{{\rm{co}}{{\rm{m}}_l}})}^T}\mathit{\boldsymbol{w}}_{\rm{b}}^{{\rm{co}}{{\rm{m}}_l}}/[(N + M) - 4]} $	（11）

辨识统计量越小，代表该组合下的伪距与正常工作通道的伪距间的一致性越好^[29]，其组合更可能对应真实信号估计伪距，进而可确定各个通道的真实信号成分。即

$ {\mathit{\boldsymbol{\rho }}^{{\rm{co}}{{\rm{m}}_l}}} = \left\{ {\begin{array}{*{20}{l}} {{\rm{ 全真实伪距组合 }}}&{T_{{\rm{dis}}}^{{\rm{co}}{{\rm{m}}_l}} = T_{{\rm{dis}}}^{\min }}\\ {{\rm{ 非全真实伪距组合 }}}&{T_{{\rm{dis}}}^{{\rm{co}}{{\rm{m}}_l}} > T_{{\rm{dis}}}^{\min }} \end{array}} \right. $	（12）

$ T_{{\rm{dis}}}^{{\rm{co}}{{\rm{m}}_l}} = T_{{\rm{dis}}}^{\min } = {\rm{min}}(T_{{\rm{dis}}}^{{\rm{co}}{{\rm{m}}_{\rm{1}}}},T_{{\rm{dis}}}^{{\rm{co}}{{\rm{m}}_{\rm{2}}}}, \cdots ,T_{{\rm{dis}}}^{{\rm{co}}{{\rm{m}}_{\rm{2}}}N}) $	（13）

利用辨识出的真实伪距组合ρ^com_l重新定位，即可获得接收机的真实位置P^cor与接收机钟差τ_rec，即

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{X}}^{{\rm{cor}}}} = \left[ {\begin{array}{*{20}{l}} {{P^{{\rm{cor}}}}}\\ {{\tau _{{\rm{rec}}}}} \end{array}} \right] = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{arg}}\mathop {{\rm{min}}}\limits_{{\mathit{\boldsymbol{X}}^{{\rm{cor}}}}} \sum\limits_{i = 1}^{N + M} {(\left\| {\mathit{\boldsymbol{P}}_i^{{\rm{sat}}} - {\mathit{\boldsymbol{P}}^{{\rm{cor}}}}} \right\| - \rho _i^{{\rm{co}}{{\rm{m}}_l}} + {\tau _{{\rm{rec}}}} - \tau _i^{{\rm{sat}}})} } \end{array} $

（14）

式中：P_i^sat为第i通道对应卫星的几何位置；ρ_i^com_l为对应的伪距；τ_i^sat为对应卫星的时钟误差。该目标泛函取最小值时，则对应接收机的估计真实导航解算结果X^cor。

具体的辨识算法流程图如图 4所示。对于第l个伪距组合，辨识统计量T_dis^com_l的获取依赖于前一历元对应伪距组合的定位结果，并以此为基础计算关系矩阵与伪距观测残差余量，进而获得对应的伪距残差w_b^com_l与T_dis^com_l。对比2^N个伪距组合对应的辨识量后，最终利用选取的真实信号伪距组合修正接收机的定位解算结果。需要注意的是，由于伪距组合数目随受欺骗影响卫星数增加而快速增加，辨识算法需要计算2^N次统计量。因此在具体实现时，可以充分利用每个通道估计的2个伪距必然为一真一假这一条件，在获得至少5个真实伪距后，每次新增加一个通道的一个伪距，如果检测量显著增大，则该新增伪距对应欺骗信号，另一伪距则为真实伪距，反之，则新增伪距为真实伪距，另一伪距为欺骗伪距，优化后辨识统计量的计算量为2⁵+(N-5)。

4 仿真验证与分析

为了测试本文所提方法的工作性能，利用欺骗模拟数据对所提方法进行测试。欺骗数据利用欺骗模拟仿真器在实测真实卫星信号的基础上模拟产生^[30]。其中，实测数据在南京航空航天大学自动化学院1号楼顶，利用东方联星NS210M接收前端采集GPS L1信号并存储在计算机中，采集系统如图 5所示，主要参数如表 2所示。

为测试所提方法中欺骗辨识算法的有效性，选取接收卫星信号较少情况时的实测数据生成欺骗仿真数据，以模拟较恶劣的欺骗环境。实测数据选取仅包含5颗可见星信号情况，如图 6所示，PRN2、6、9、17和19卫星可见。

为测试所提方法对欺骗信号与真实信号参数的估计能力，欺骗仿真采用较难以处理的牵引式欺骗，对其中任意2颗卫星信号进行欺骗攻击，具体地，欺骗仿真器在第10000 ms对PRN6与PRN17卫星进行牵引式欺骗攻击，牵引速率(相对码延时率)分别为24 m/s与40 m/s，欺骗信号与真实信号幅值比R^ar=2(功率比J/S=6.02 dB)。其他可见星PRN2、9和19卫星信号则未受欺骗干扰影响。

所生成的欺骗仿真信号分别用软件定义的传统接收机与论文所提抗欺骗方法进行解算。传统接收机跟踪环路参数如表 2所示；基于EKF的真实与欺骗信号参数估计环路在检测到欺骗信号后开始工作，假设欺骗检测算法检测时间为500 ms，估计环路在欺骗干扰后500 ms启动。

图 7与图 8分别展示了基于EKF的真实与欺骗信号参数估计环路对欺骗干扰下PRN6与PRN17卫星的参数估计结果。对于PRN6卫星信号，欺骗信号于第10000 ms引入接收机，且其牵引速率24 m/s≈0.0819 chip/s，对应的预设相对欺骗码延时α=α^sp+α^au在第10500 ms时应约为0.041 chip，第15000 ms时达到0.409 chip。图 7(a)中，估计环路在10500 ms开始工作，估计α曲线基本与预设欺骗码延时变化吻合，增长率约0.0793 chip/s。图 7(b)中，A^au与A^sp的平均估计结果分别为2915与5697，A^sp/A^au的平均估计结果为1.956，与预设幅值比2近似，误差0.044。类似地，参数估计环路对PRN17卫星信号的估计结果也符合欺骗数据的设置。图 8中预设欺骗码延时在10500 ms时为0.068 chip，在15000 ms时达到0.682 chip，环路估计的α与该预设α变化一致，真实信号自相关最大幅值的平均估计值约3310，而欺骗信号的自相关最大幅值约5997，A^sp/A^au的平均估计结果为1.813，与预设值相差0.187。整体上，所提估计方法能够比较准确地估计出真实信号与欺骗信号的幅值A^au、A^sp与码延时α^sp、α^au参数。

各欺骗通道在真实与欺骗信号参数估计环路估计信息的辅助下，得到2个估计伪距，对应共4种伪距组合：

$ \left\{ {\begin{array}{*{20}{l}} {{\rho ^{{\rm{co}}{{\rm{m}}_{\rm{1}}}}}:\{ \rho _{{\rm{PR}}{{\rm{N}}_{{\rm{17}}}}}^1,\rho _{{\rm{PR}}{{\rm{N}}_{\rm{6}}}}^1,{\rho _{{\rm{PR}}{{\rm{N}}_{{\rm{19}}}}}},{\rho _{{\rm{PR}}{{\rm{N}}_{\rm{9}}}}},{\rho _{{\rm{PR}}{{\rm{N}}_{\rm{2}}}}}\} }\\ {{\rho ^{{\rm{co}}{{\rm{m}}_{\rm{2}}}}}:\{ \rho _{{\rm{PR}}{{\rm{N}}_{{\rm{17}}}}}^1,\rho _{{\rm{PR}}{{\rm{N}}_{\rm{6}}}}^2,{\rho _{{\rm{PR}}{{\rm{N}}_{{\rm{19}}}}}},{\rho _{{\rm{PR}}{{\rm{N}}_{\rm{9}}}}},{\rho _{{\rm{PR}}{{\rm{N}}_{\rm{2}}}}}\} }\\ {{\rho ^{{\rm{co}}{{\rm{m}}_{\rm{3}}}}}:\{ \rho _{{\rm{PR}}{{\rm{N}}_{{\rm{17}}}}}^2,\rho _{{\rm{PR}}{{\rm{N}}_{\rm{6}}}}^1,{\rho _{{\rm{PR}}{{\rm{N}}_{{\rm{19}}}}}},{\rho _{{\rm{PR}}{{\rm{N}}_{\rm{9}}}}},{\rho _{{\rm{PR}}{{\rm{N}}_{\rm{2}}}}}\} }\\ {{\rho ^{{\rm{co}}{{\rm{m}}_{\rm{4}}}}}:\{ \rho _{{\rm{PR}}{{\rm{N}}_{{\rm{17}}}}}^2,\rho _{{\rm{PR}}{{\rm{N}}_{\rm{6}}}}^2,{\rho _{{\rm{PR}}{{\rm{N}}_{{\rm{19}}}}}},{\rho _{{\rm{PR}}{{\rm{N}}_{\rm{9}}}}},{\rho _{{\rm{PR}}{{\rm{N}}_{\rm{2}}}}}\} } \end{array}} \right. $

对应的各伪距组合下的辨识检测量分别为T_dis^com₁、T_dis^com₂、T_dis^com₃、T_dis^com₄。辨识检测量曲线如图 9所示：在第10500~30000 ms之间，检测量T_dis^com₁相比其他3组检测量的取值始终最小，均值为4.05 m。该检测统计量所对应的伪距组合为：{ρ_PRN17¹, ρ_PRN6¹, ρ_PRN19, ρ_PRN9, ρ_PRN2}，与预设欺骗场景一致。利用参数估计环路辅助辨识算法可正确辨识出该全真实伪距组。另外，随着欺骗时间增加，欺骗信号的码延时相位偏离真实信号更远，对应地，欺骗伪距与真实伪距的不一致性也越加严重，这体现在了包括欺骗伪距组合的辨识检测量T_dis^com₂、T_dis^com₃、T_dis^com₄的变化情况：随时间逐渐增加，远远大于全真实伪距组合下的检测量T_dis^com₁。例如，在前2000 ms，T_dis^com₂取值均值为114.76 m，到后2000 ms(28000~30000 ms)，T_dis^com₂增加到约280 m，远高于真实伪距组检测量T_dis^com₁。

图 10为传统接收机与抗欺骗接收机在东北天坐标系下的定位偏差曲线，取接收机未受欺骗干扰阶段(0~10000 ms)定位均值为坐标原点。

传统接收机在欺骗攻击下解算坐标迅速偏离，受可见星分布与受欺骗卫星分布的影响，欺骗导致的定位偏差主要体现在天向与北向定位结果上，在第35000 ms时，天向与北向位置偏差值分别为2172.16 m与4403.48 m，三维位置偏差高达4919.35 m。而抗欺骗接收机于15000 ms时完成欺骗估计与辨识过程，并迅速修正定位结果，使得接收机回归正常位置。在第35000 ms时，抗欺骗接收机的定位偏差仅9.27 m，与未受欺骗干扰时定位结果基本一致。

欺骗参数估计算法是所提抗欺骗方法的核心算法，为测试参数估计算法性能，测试中分别调整欺骗信号的相对码延时α=0.3, 0.5, 0.7, 0.8, 0.9, 1.0 chip与幅值比R^ar=1, 2, 3, 4, 5, 6，牵引速率0 m/s(退化为非牵引式欺骗)，观察欺骗估计算法对码延时和幅值比的估计误差，如图 11所示。

图 11反映了对码延时估计误差的变化趋势，随相对欺骗码延时的增大，幅值比估计误差与码延时估计误差均增加。在相对欺骗码延时0.3~0.7 chip间，幅值与码延时估计精度均较高，幅值比估计误差小于0.2，相对码延时估计误差小于0.1 chip，考虑幅值比估计误差大小不影响伪距估计结果，也不影响定位结果，从图 11(b)中可见，在预设相对欺骗码延时介于0.3~0.9 chip之间，真实与欺骗信号参数估计环路关于码延时估计误差的RMSE(Root Mean Square Error)低于0.1 chip。

图 12(a)与图 12(b)所示为预设相对欺骗码延时为0.3、0.5、0.7 chip时，幅值比估计误差与相对码延时估计误差随预设欺骗/真实信号幅值变化曲线。随预设幅值比从1增加到5，幅值估计与码延时估计具有较高估计精度，码延时估计误差低于0.1 chip，幅值比估计误差低于1.65。

简言之，对牵引式欺骗以及非牵引式欺骗，当欺骗与真实信号相对码延时为0.3~0.9 chip，幅值比在1~5间，所提参数估计算法能够利用九路相关器实现较高精度的欺骗信号与真实信号参数估计，幅值比估计误差小于1.65，相对码延时估计误差小于0.1 chip。

5 结论

本文研究了一种基于卡尔曼滤波技术辅助RAIM算法的抗欺骗技术，其具有仅利用九路相关器估计混合信号中真实与欺骗信号的幅值与码延时参数的能力，并能够辨识出真实信号。所提方法使得接收机在预设的多星欺骗场景下，有效修正定位偏差至10 m内；同时对于预设的相对欺骗码延时介于0.3~0.9 chip的欺骗信号，跟踪环路对码延时估计的误差一般低于0.1 chip.

未来工作集中将针对估计算法中相关器的选择进行更深入的研究与分析；利用其他导航信息辅助跟踪环路与辨识算法，改进参数估计精度，提高辨识算法的效率。