CN103926596A

CN103926596A - 一种基于粒子滤波的稳健gnss抗欺骗方法

Info

Publication number: CN103926596A
Application number: CN201410169614.2A
Authority: CN
Inventors: 孟维晓; 巩紫君; 韩帅; 罗德巳
Original assignee: Harbin Institute of Technology
Current assignee: Harbin Institute of Technology
Priority date: 2014-04-25
Filing date: 2014-04-25
Publication date: 2014-07-16
Anticipated expiration: 2034-04-25
Also published as: CN103926596B

Abstract

一种基于粒子滤波的稳健GNSS抗欺骗方法，本发明涉及GNSS抗欺骗方法。是要解决欺骗卫星产生偏差，无法抵抗转发式欺骗以及抗欺骗无法应用于已有的系统，产生成本高而提出的GNSS抗欺骗方法，该方法是通过1、得到1时刻状态向量粒子；2、更新k时刻状态向量粒子；3、得到观测伪距4、更新观测向量5、计算出未归一化的粒子权重6、计算观测向量的模值的最小值θ；7、根据修正8、计算k时刻定位结果9、得到剩余的有效向量粒子的数目M_eff；10、计算出11、将带入到到步骤二等步骤实现的。本发明应用于GNSS抗欺骗领域。

Description

一种基于粒子滤波的稳健GNSS抗欺骗方法

技术领域

本发明涉及一种基于粒子滤波的稳健GNSS抗欺骗方法。

背景技术

随着LBS(Location Based Service，基于位置的服务)需求的不断上升和接收机价格的不断下降，GNSS(Global Navigation Satellite System，全球导航卫星系统)已经广泛应用于军用和民用市场。GNSS不仅能够提供精确的定位和授时服务，而且被普遍应用与很多新兴的无线系统，比如智能电网和CDMA2000通信系统。但随着接受度的增高，GNSS的弱点也严重暴露了出来，尤其是对欺骗式攻击，几乎毫无抵抗力。欺骗式攻击就是利用一个无线电设备，模拟产生卫星的欺骗式导航信号。欺骗提供错误的伪距信息或卫星位置信息，从而误导接收机，使它产生偏差很大的定位结果。

十几年前，很少有人研究怎样抵抗欺骗式攻击的问题，因为那时候，这种攻击的实现过于复杂。但随着集成电路技术的快速发展，欺骗式攻击已经具备了条件，所以必须被纳入学者的考虑范围之内。

最早的抗欺骗方法应用在美国军事上，利用GPS军用信号的加密伪随机码，只要干扰方不知道这个扩频码，就无法模拟欺骗信号。但这一方法无法估计到民用用户，而且无法抵抗转发式欺骗。Logan Scott提出了一种新的信号架构来实现抗欺骗，K.Wesson也讨论了这一问题，但这都需要改变信号体制，无法应用于已有的系统。C.J.Wullems和X.Jiang也提出了自己的方法，但他们的问题是需要额外的硬件设备，提高了成本。

发明内容

本发明的目的是为了解决欺骗卫星提供错误的伪距信息误导接收机，使它产生偏差很大的定位结果，无法抵抗转发式欺骗以及抗欺骗需要改变信号体制，无法应用于已有的系统，需要额外的硬件设备，提高了成本的问题而提出的一种基于粒子滤波的稳健GNSS抗欺骗方法。

上述的发明目的是通过以下技术方案实现的：

步骤一、确定用户的初始状态，在用户初始状态的领域内均匀分布m个初始状态向量粒子同时设置迭代次数的计数值为k＝1；根据用户初始状态向量粒子可以得到k＝1时刻用户状态向量粒子

x_{1}^{(m)} = x_{0}^{(m)} + v_{1}^{(m)} + f_{1}^{(m)},

其中，

v_{1} = {[{\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, {\overset{\cdot}{δ} t}_{U, 1}]}^{T}, {\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, \overset{\cdot}{δ} t_{U, 1}

分别表示x_U,1,y_U,1,z_U,1,δt_U,1在单位时隙中的变化量；表示第k＝1时刻的状态向量粒子的过程噪声，服从均值为零的高斯分布；x_U,1,y_U,1,z_U,1代表用户在k＝1个时刻状态向量粒子的三维位置,δt_U,1用户在k＝1个时刻状态与第n颗卫星的时钟偏差；

步骤二、当k≠1时，根据k-1时刻的用户状态向量粒子更新m个在第k时刻状态向量粒子

x_{k}^{(m)} = x_{k - 1}^{(m)} + v_{k}^{(m)} + f_{k}^{(m)}; k = 2.3 . . . . .

其中，状态向量粒子(x_U,k,y_U,k,z_U,k)代表用户在k个时刻状态向量粒子的三维位置；δt_U,k用户在k个时刻状态与第n颗卫星的时钟偏差；模型输入为：

v_{k}^{(m)} = {[{\overset{\cdot}{x}}_{U, k}, {\overset{\cdot}{y}}_{U, k}, {\overset{\cdot}{z}}_{U, k}, {\overset{\cdot}{δ} t}_{U, k}]}^{T}, {\overset{\cdot}{x}}_{U, k}, {\overset{\cdot}{y}}_{U, k}, {\overset{\cdot}{z}}_{U, k}, \overset{\cdot}{δ} t_{U, k}

分别表示x_U,k,y_U,k,z_U,k,δt_U,k在单位时隙中的变化量；表示用户在第k时刻的状态向量粒子的过程噪声，服从均值为零的高斯分布；

步骤三、根据步骤二计算的计算出用户与第n颗卫星之间的观测伪距由下式给出：

是接收机和卫星之间的观测伪距，(x_n,k,y_n,k,z_n,k)(n＝1,2,...,N)代表第n颗卫星在第k个时刻的位置坐标，c表示光速；

步骤四、根据步骤三得到的观测伪距更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T},

其中，用户与第n颗卫星的伪距观测误差为ρ_n,k表示在第k个时刻，第n颗卫星与用户之间的真实伪距，N为可见卫星的个数；

步骤五、根据更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T},

计算未归一化的粒子权重

w_{k}^{(m)} = \frac{1}{M} {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}]

其中，M为在第k时刻的状态向量粒子数，N为可见卫星的个数，σ为的标准差；

步骤六、根据未归一化的粒子权重计算观测向量的模值的最小值θ：

\begin{matrix} θ = - 2 σ^{2} (\ln (M \max (w_{k}^{(m)})) + \frac{N}{2} \ln (2 π σ^{2})) \\ = \min {| | y_{k}^{(m)} | |}^{2} \end{matrix};

步骤七、判断θ与γ的大小，如果θ>γ，根据修正得到：

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix}

其中，b为误差下界，d为误差上界，为观测误差，P_f(·)为误差抑制函数；

然后根据修正后的再重新计算未归一化的粒子权重如果θ≤γ，则直接跳到步骤八；

步骤八、根据未归一化的粒子权重计算k时刻定位结果

{\overset{&OverBar;}{w}}_{k}^{(m)} = \frac{w_{k}^{(m)}}{Σ_{m = 1}^{M} w_{k}^{(m)}},

{\hat{x}}_{k}^{MMS} = Σ_{m = 1}^{M} {\overset{&OverBar;}{w}}_{k}^{(m)} x_{k}^{(m)};

其中，为归一化粒子权重；

步骤九、将步骤二中共得到M个粒子在k时刻的状态向量粒子根据步骤五中的未归一化的粒子权重计算得到剩余的有效向量粒子的数目M_eff：

M_{eff} = \frac{1}{Σ_{m = 1}^{M} {(w_{k}^{(m)})}^{2}},

如果则转到步骤十一，如果跳到步骤十，进行重采样；

步骤十、如果采取M个在k时刻的状态向量粒子对集合重采样，采到的概率为归一化粒子权重

步骤十一、如果利用步骤二中的概率为的在k时刻状态向量粒子将k时刻加1后得到用户在k+1时刻的用户状态向量粒子跳到步骤二，开始下一个时刻的定位过程；即完成了一种基于粒子滤波的稳健GNSS抗欺骗方法。

发明效果

本发明提出了一种基于粒子滤波的抗欺骗方法，能够在不增加任何硬件资源、不改变信号体制的条件下实现欺骗信号的侦查和消除。从而确保了定位精度的准确，保护了国家和人民的利益。该发明可以有效侦测和抑制欺骗卫星的干扰，图2可以看到，可见卫星越多、欺骗卫星给出的伪距误差越大，则检测效果越好。在误差达到40m时，检测概率可以达到70％以上。为了证明该算法对欺骗卫星干扰的抑制，下面对传统的最小二乘法、传统的粒子滤波算法和本文的算法做了比较。图5、图6可以看到，仿真时间为250秒，前120秒没有干扰，之后加入一颗欺骗卫星，仿真时间为250秒0秒到120秒没有干扰，三种方法的定位误差差不多，在第120秒加入干扰后，传统的最小二乘法和粒子滤波方法都产生了很大的定位误差，但本文的算法却还是保持了很好的定位精度。转发式的干扰，其主要机理是提供假的伪距，从而欺骗接收机。而本发明的本质，就是检测伪距的残差，通过伪距的异常来判断是否存在欺骗卫星。可见，本发明可以检测一切提供假伪距的欺骗方法，转发式干扰自然也包括在内。本发明不需要新的信号架构，即不需要改变信号体制，可以直接应用于已有的系统。这一点从后面的分析过程中可以看出，该方法使用的参数都是现有系统可以提供的

本发明直接应用于已有系统，不需要任何额外的硬件，只需要修改解算过程，即对接收机的软件做调整即可。而现在多数其他的抗欺骗式干扰的方法是需要额外的硬件的，所以本发明相对来说成本低。

附图说明

图1是具体实施方式一提出的一种基于粒子滤波的稳健GNSS抗欺骗方法流程图；

图2是具体实施方式一提出的发现欺骗卫星的概率P_D与卫星数N的关系图，其中横坐标△P_f为欺骗卫星提供的伪距和卫星真实伪距的偏差，纵坐标为发现欺骗卫星的概率P_D；为自由度为N＝5时的发现欺骗卫星的概率P_D，为自由度为N＝8时的发现欺骗卫星的概率P_D；

图3是具体实施方式四提出的自由度N＝5时g(θ)与g_f(θ)的关系图；其中，横坐标为监测变量θ，纵坐标PDF of θ为监测变量θ的概率密度函数值；为不存在欺骗卫星，监测变量θ的概率密度函数g(θ)值变化曲线；为存在欺骗卫星，则监测变量θ的概率密度函数g_f(θ)值的变化曲线；

图4是具体实施方式四提出的N＝8时g(θ)与g_f(θ)的关系图；其中，横坐标为监测变量θ，纵坐标PDF of θ为监测变量θ的概率密度函数值；为不存在欺骗卫星，监测变量θ的概率密度函数g(θ)值变化曲线；为存在欺骗卫星，则监测变量θ的概率密度函数g_f(θ)值的变化曲线；

图5是具体实施方式一提出的N＝5时不同算法的定位误差与定位时间曲线图；其中为传统的最小二乘法，为传统的抗干扰方法，为本发明提出的抗干扰算法；

图6是具体实施方式一提出的N＝8时不同算法的定位误差与定位时间曲线图；其中为传统的最小二乘法，为传统的抗干扰方法，为本发明提出的抗干扰算法。

具体实施方式

具体实施方式一：本实施方式的一种基于粒子滤波的稳健GNSS抗欺骗方法，具体是按照以下步骤制备的：

步骤一、完成两个变量的初始化，首先将确定用户的初始状态，在用户初始状态的领域内均匀分布m个初始状态向量粒子均匀分布在x₀领域内；同时设置迭代次数的计数值为k＝1；根据用户初始状态向量粒子可以得到k＝1时刻用户状态向量粒子

x_{1}^{(m)} = x_{0}^{(m)} + v_{1}^{(m)} + f_{1}^{(m)},

其中，

v_{1} = {[{\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, {\overset{\cdot}{δ} t}_{U, 1}]}^{T}, {\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, \overset{\cdot}{δ} t_{U, 1}

分别表示x_U,1,y_U,1,z_U,1,δt_U,1在单位时隙中的变化量；表示第k＝1时刻的状态向量粒子的过程噪声，服从均值为零的高斯分布；x_U,1,y_U,1,z_U,1代表用户在k＝1个时刻状态粒子的三维位置,δt_U,1用户在k＝1个时刻状态与第n颗卫星的时钟偏差；

x_{k}^{(m)} = x_{k - 1}^{(m)} + v_{k}^{(m)} + f_{k}^{(m)}; k = 2,3 . . . . .

是粒子滤波的动态处理模型，以适应GNSS接收机的运动；其中，状态向量粒子也就是用户在k时隙的状态，也是最终的定位结果；所有坐标都是在ECEF(Earth-Centered and Earth-Fixed，地心固定坐标系)坐标系下给出，(x_U,k,y_U,k,z_U,k)代表用户在k个时刻状态粒子的三维位置；δt_U,k用户在k个时刻状态与第n颗卫星的时钟偏差；模型输入为：分别表示x_U,k,y_U,k,z_U,k,δt_U,k在单位时隙中的变化量；表示用户在第k时刻的状态向量粒子的过程噪声，服从均值为零的高斯分布；

步骤四、根据步骤三得到的观测伪距更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T}

步骤五、根据更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T}

计算未归一化的粒子权重

w_{k}^{(m)} = \frac{1}{M} {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}]

\begin{matrix} θ = - 2 σ^{2} (\ln (M \max (w_{k}^{(m)})) + \frac{N}{2} \ln (2 π σ^{2})) \\ = \min {| | y_{k}^{(m)} | |}^{2} \end{matrix}

步骤七、判断θ与γ的大小，如果θ>γ，根据修正得到：

其中

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix}

然后根据修正后的再重新计算未归一化的粒子权重使得θ≤γ；

如果θ≤γ，则直接跳到步骤八；

步骤八、根据未归一化的粒子权重计算k时刻定位结果

{\overset{&OverBar;}{w}}_{k}^{(m)} = \frac{w_{k}^{(m)}}{Σ_{m = 1}^{M} w_{k}^{(m)}},

{\hat{x}}_{k}^{MMS} = Σ_{m = 1}^{M} {\overset{&OverBar;}{w}}_{k}^{(m)} x_{k}^{(m)};

其中，为归一化粒子权重；

步骤九、将步骤二中共得到M个粒子在k时刻的状态向量粒子由于存在粒子退化现象，部分状态粒子将失效，剩余的有效状态粒子的数目M_eff，根据步骤五中的未归一化的粒子权重计算得到剩余的有效向量粒子的数目M_eff：

M_{eff} = \frac{1}{Σ_{m = 1}^{M} {(w_{k}^{(m)})}^{2}},

如果则转到步骤十一，如果跳到步骤十，进行重采样；

步骤十、如果采取M个在k时刻的状态向量粒子对集合重采样，采到的概率为归一化粒子权重重新采到的归一化粒子权重直至使得为止；

步骤十一、如果利用步骤二中的概率为的在k时刻状态向量粒子将k时刻加1后得到用户在k+1时刻的用户状态向量粒子跳到步骤二，开始下一个时刻的定位过程；如图1即完成了一种基于粒子滤波的稳健GNSS抗欺骗方法。

本实施方式效果：

本实施方式提出了一种基于粒子滤波的抗欺骗方法，能够在不增加任何硬件资源、不改变信号体制的条件下实现欺骗信号的侦查和消除。从而确保了定位精度的准确，保护了国家和人民的利益。本实施方式可以有效侦测和抑制欺骗卫星的干扰，图2可以看到，可见卫星越多、欺骗卫星给出的伪距误差越大，则检测效果越好。在误差达到40m时，检测概率可以达到70％以上。为了证明该算法对欺骗卫星干扰的抑制，下面对传统的最小二乘法、传统的粒子滤波算法和本文的算法做了比较。图5、图6可以看到，仿真时间为250秒，前120秒没有干扰，之后加入一颗欺骗卫星，仿真时间为250秒0秒到120秒没有干扰，三种方法的定位误差差不多，在第120秒加入干扰后，传统的最小二乘法和粒子滤波方法都产生了很大的定位误差，但本文的算法却还是保持了很好的定位精度。转发式的干扰，其主要机理是提供假的伪距，从而欺骗接收机。而本实施方式的本质，就是检测伪距的残差，通过伪距的异常来判断是否存在欺骗卫星。可见，本实施方式可以检测一切提供假伪距的欺骗方法，转发式干扰自然也包括在内。本实施方式不需要新的信号架构，即不需要改变信号体制，可以直接应用于已有的系统。这一点从后面的分析过程中可以看出，该方法使用的参数都是现有系统可以提供的

本实施方式直接应用于已有系统，不需要任何额外的硬件，只需要修改解算过程，即对接收机的软件做调整即可。而现在多数其他的抗欺骗式干扰的方法是需要额外的硬件的，所以本发明相对来说成本低。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤三中根据步骤二计算的计算出用户与第n颗卫星之间的观测伪距由下式给出：

(1)通过用户在第k时刻的状态向量粒子计算得到观测伪距

其中，(x_n,k,y_n,k,z_n,k)(n＝1,2,...,N)代表在第n颗卫星在k个时刻的位置坐标；(x_U,k,y_U,k,z_U,k,δt_U,k)代表用户在第k时刻的状态向量粒子的状态，包括用户在第k时刻的状态向量粒子三维位置(x_U,k,y_U,k,z_U,k)和相对于第n颗卫星的时钟偏差δt_U,k；

(2)真实伪距ρ_n,k由接收机的伪距测量电路直接测量得到，假设可见卫星的数目为N，信号捕获和跟踪的影响被忽略；在第k个时刻，测得用户和第n颗卫星之间的观测伪距与真实伪距ρ_n,k的关系为：

ε_n,k表示在k时刻接收机测量与第n颗卫星的伪距时产生的测量误差。其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤五中根据更新观测向量计算未归一化的粒子权重的过程为：

(1)建立离散时间的粒子滤波模型为：

\{\begin{matrix} {x_{k}}^{(m)} = {x^{(m)}}_{k - 1} + {v^{(m)}}_{k} + {f_{k}}^{(m)} \\ {y_{k}}^{(m)} = h^{(m)} ({x_{k}}^{(m)}) + {ϵ^{(m)}}_{k} \end{matrix}

(由y_k ^(m)求出由判断是否需要重采样，如果重采样，也就是重新计算在步骤十中体现)

观测向量为代表了测量得到的伪距与理论计算得到的伪距之间的差值；ε^(m) _k表示第k时刻的状态向量粒子的观测噪声；h^(m)(·)代表与y_k ^(m)之间的函数关系；

(2)计算未归一化的粒子权重的过程为：

其中，未归一化的粒子权重可由下式得到：

w_{k}^{(m)} = \frac{p (y_{k}^{(m)} | x_{k}^{(m)}) \hat{p} (x_{k}^{(m)} | x_{k - 1}^{(m)})}{q (x_{k}^{(m)} | x_{k - 1}^{(m)}, y_{k}^{(m)})} w_{k - 1}^{(m)}

上式中，q(·)为重要密度函数(这个是为了推导方便，人为定义的一个函数，没有明确的意义)，每一个状态向量粒子都来自对q(·)的完备采样，为已知时的概率密度函数，为已知时的概率密度函数；为了得到方差最小的未归一化粒子权重，应该选取的后验概率密度即已知和时的概率密度，但在大多数实际情况中是无法实现的，可以采用近似表示，得到了未归一化粒子权重的表达式：

w_{k}^{(m)} = p (y_{k}^{(m)} | x_{k}^{(m)}) w_{k - 1}^{(m)}

其中，经过重采样后，我们可以认为

w_{k - 1}^{(m)} = \frac{1}{M}

p (y_{k}^{(m)} | x_{k}^{(m)}) = {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}]

不失一般性地假设第N颗卫星是欺骗卫星，从而得到了未归一化粒子权重的最终表达：

\begin{matrix} w_{k}^{(m)} = p (y_{k}^{(m)} | x_{k}^{(m)}) w_{k - 1}^{(m)} \\ = \frac{1}{M} {(\frac{1}{\sqrt{2 π σ^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}] \end{matrix}

归一化后为：

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤六中计算观测向量的模值的最小值θ：

假设有多于四颗可见的真实卫星，还有一颗欺骗卫星；这颗欺骗卫星会给接收机提供一个自身坐标(x_f,k,y_f,k,z_f,k)和假的伪距从而提出了对欺骗的侦测和抑制方法，由下式给出：

{\hat{ρ}}_{f, k} = ρ_{f, k} + {Δρ}_{f, k}

ρ_f,k是接收机和卫星之间的真实伪距，是欺骗卫星提供的假伪距，△ρ_f,k是二者的偏差；一般而言，定位误差与|△ρ_f,k|正相关；(本公式是理论上的推导，说明本方法的合理性)

前面已经提到，就是测量伪距与理论计算得到的伪距之间的差值；在没有欺骗卫星的条件下，该向量中的每一个数都不会很大；如果存在欺骗卫星，对应的伪距误差就会很大，从而导致的二范数很大，也就导致很大；所以，通过监测来侦测欺骗卫星，具体如下：

\max (w_{k}^{(m)}) = \frac{1}{M} {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{\min (| | y_{k}^{(m)} | |)}{2 σ^{2}}]

根据定义监测变量θ为：

\begin{matrix} θ = - 2 σ^{2} (\ln (M \max (w_{k}^{(m)})) + \frac{N}{2} \ln (2 π σ^{2})) \\ = \min {| | y_{k}^{(m)} | |}^{2} \end{matrix};

当存在欺骗卫星时且自由度N＝5时，θ服从非中心卡方分布，其概率密度函数为：

g_{f} (θ) = \{\begin{matrix} \frac{1}{2 \sqrt{2 π σ^{2} θ}} \exp (- \frac{θ + λ}{2 σ^{2}}) \cdot (\exp (\frac{\sqrt{θλ}}{σ^{2}}) + \exp (- \frac{\sqrt{θλ}}{σ^{2}})), & θ &GreaterEqual; 0 \\ 0, & θ &GreaterEqual; 0 \end{matrix}

在N>5的条件下，g_f(θ)可以近似为：

g_{f} (θ) = \{\begin{matrix} \frac{1}{2 σ^{2}} {(\frac{θ}{λ})}^{\frac{N - 6}{4}} \exp (- \frac{θ + λ}{2 σ^{2}}) I_{\frac{N - 6}{2}} (\frac{\sqrt{θλ}}{σ^{2}}), & θ &GreaterEqual; 0 \\ 0, & θ &GreaterEqual; 0 \end{matrix}

I_T(·)为T阶修正的贝塞尔函数，

如果不存在欺骗卫星，则θ服从中心卡方分布，自由度为N-4：

g (θ) = \{\begin{matrix} \frac{1}{2 σ^{2} Γ (\frac{N}{2})} {(- \frac{θ}{{2 σ}^{2}})}^{\frac{N - 2}{2}} \exp (- \frac{θ}{2 σ^{2}}), & θ &GreaterEqual; 0 \\ 0, & θ < 0 \end{matrix}

其中，Γ(·)是gamma函数，λ为卡方分布的非中心度参数，其取值与可见卫星个数N有关，也与伪距有关；由下面的方法得出：

λ = {(b_{k})}^{T} (I - G_{k} {(G_{k}^{T} G_{k})}^{- 1} G_{k}^{T}) b_{k}

b_k＝[0 0 ... b_s,k ... 0]^T

G_{k} = [\begin{matrix} 1_{1, k} & 1 \\ 1_{2, k} & 1 \\ . & . \\ . & . \\ . & . \\ 1_{N, k} & 1 \end{matrix}]

在第k个时刻，这里1_n,k是当前接收机指向第n颗卫星的单位方向矢量，I表示单位矩阵，b_s,k为第n颗卫星的附加伪距，b_k为对应的附加伪距矩阵，G_k表示接收机与所有可见卫星的单位方向矩阵，T代表矩阵的转置，是G_k的转置；

可见，存在欺骗卫星的情况下和不存在欺骗卫星的情况下，θ服从不同的分布，可以通过该分布来判断是否存在欺骗卫星；因此，可以设置一个门限γ，当θ>γ时认为存在欺骗卫星，反之则认为没有；依据T阶修正的贝塞尔函数I_T(·)说明本方法的合理性如图3、4所示。其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤七中判断θ与γ的大小，如果θ>γ，根据修正得到

其中

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix},

P_f(·)为误差抑制函数，b为误差下界，d为误差上界，为观测误差；

然后根据修正后的再重新计算未归一化粒子权重

当不存在欺骗卫星时，该方法有一定的概率发生误判，这一误判的概率被称为虚警概率P_f，与门限γ有关；虚警概率P_f与门限γ的关系由下面的式子给出：

P_{f} = {&Integral;}_{γ}^{+ \infty} g (θ) dθ,

相应的，发现概率为：

P_{D} = {&Integral;}_{γ}^{+ \infty} g_{f} (θ) dθ;

其中，g(θ)：不存在欺骗卫星的情况下，监测变量θ的概率密度函数；P_D：发现欺骗卫星的概率；g_f(θ)：存在欺骗卫星的情况下，监测变量θ的概率密度函数；

在侦测到欺骗卫星后，可以通过基于稳健估计的改进方法来对其进行抑制；稳健估计的目的在于，在一些普遍假设的条件下，给出最贴近真实值的估计；传统的估计方法，比如最小二乘法，认为每一个观测数据都服从先验条件，所以，某一个错误的观测值会带来很大的误差；而稳健估计则不同，它会去除某些异常观测值，从而消除其影响；最大似然类型的估计(M-estimate)，高阶统计的线性组合(L-estimate)，等，是基本的稳健估计方法；其中，最大似然比估计总是被用在参数估计上；

本发明基于传统的稳健估计方法，提出了改进，具体处理方法如下：

首先定义误差抑制函数P_f(·)：

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix}

其中，b为误差下界，d为误差上界，e为观测误差，在计算得到未归一化粒子权重后，对观测误差进行修正，得到方法如下：

经过这一处理后：

(1)若误差值的绝对值超过d，则说明对应的卫星是一颗假卫星，将它修正为0，即它就不会对后续计算产生任何影响了，也就是说不会影响到定位结果，所以假卫星被完全抑制，

(2)若的绝对值在b和d之间，则修正后减小为b，即

(3)若的绝对值小于b，则其大小则保持不变，即然后再重新计算未归一化粒子权重，得到归一化粒子权重这一部分的作用在于，抑制了欺骗卫星对定位结果的干扰。其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤八中根据未归一化的粒子权重计算k时刻定位结果推导过程为：

考虑最小均方估计，则对于的估计为

\begin{matrix} {\hat{x}}_{k}^{MMS} = Σ_{m = 1}^{M} {&Integral;}_{- \infty}^{+ \infty} y_{k}^{(m)} p (x_{k}^{(m)} | y_{k}^{(m)}) d x_{k}^{(m)} \\ = Σ_{m = 1}^{M} {\overset{&OverBar;}{w}}_{k}^{(m)} x_{k}^{(m)} \end{matrix},

即作为第k个时刻的定位结果输出。其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：步骤十中如果采取M个在k时刻的用户状态向量粒子具体过程为：

对集合中的元素进行均匀采样，依次抽取M次，得到k时刻M个新的状态向量粒子，用k时刻M个新的状态向量粒子来取代之前的k时刻的M个状态向量粒子其中M个状态向量粒子抽取的次数为M次，如果则说明有效粒子数M_eff不足，需要进行重采样，以解决粒子退化问题；所谓粒子退化问题是指，经过若干次递推后，只有少数粒子的权重很大，而大部分粒子的权重非常接近于零的现象；由于权重接近于零的粒子对状态的估计的贡献很小，所以当粒子退化严重时，有很多计算和存储资源都被浪费了。其它步骤及参数与具体实施方式一至六之一相同。

采用以下实施例验证本发明的有益效果：

实施例一

步骤一、完成两个变量的初始化，首先将确定用户的初始状态，在用户初始状态的领域内均匀分布10000个初始状态向量粒子均匀分布在x₀领域内；同时设置迭代次数的计数值为k＝1；根据用户初始状态向量粒子可以得到k＝1时刻用户状态向量粒子

x_{1}^{(m)} = x_{0}^{(m)} + v_{1}^{(m)} + f_{1}^{(m)},

其中，

v_{1} = {[{\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, {\overset{\cdot}{δ} t}_{U, 1}]}^{T}, {\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, \overset{\cdot}{δ} t_{U, 1}

x_{k}^{(m)} = x_{k - 1}^{(m)} + v_{k}^{(m)} + f_{k}^{(m)}; k = 2,3 . . . . .

是粒子滤波的动态处理模型，以适应GNSS接收机的运动；其中，状态向量粒子也就是用户在k时隙的状态，也是最终的定位结果；所有坐标都是在ECEF(Earth-Centered and Earth-Fixed，地心固定坐标系)坐标系下给出，(x_U,k,y_U,k,z_U,k)代表用户在k个时刻状态向量粒子的三维位置；δt_U,k用户在k个时刻状态与第n颗卫星的时钟偏差；模型输入为：分别表示x_U,k,y_U,k,z_U,k,δt_U,k在单位时隙中的变化量，这里的单位时隙为1秒；表示用户在第k时刻的状态向量粒子的过程噪声，服从均值为零的高斯分布；

(1)通过用户在第k时刻的状态向量粒子计算得到观测伪距

其中，是接收机和卫星之间的观测伪距，(x_n,k,y_n,k,z_n,k)(n＝1,2,...,N)代表第n颗卫星在第k个时刻的位置坐标，c表示光速；(x_U,k,y_U,k,z_U,k,δt_U,k)代表用户在第k时刻的状态向量粒子的状态，包括用户在第k时刻的状态向量粒子三维位置(x_U,k,y_U,k,z_U,k)和相对于第n颗卫星的时钟偏差δt_U,k；

ε_n,k表示在k时刻接收机测量与第n颗卫星的伪距时产生的测量误差；

步骤四、根据步骤三得到的观测伪距更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T}

其中，用户与第n颗卫星的伪距观测误差为ρn,k表示在第k个时刻，第n颗卫星与用户之间的真实伪距，N为可见卫星的个数；

步骤五、根据更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T}

计算未归一化的粒子权重

w_{k}^{(m)} = \frac{1}{M} {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}]

其中，M为在第k时刻的状态向量粒子数，取10000，N为可见卫星的个数，取5，σ为的标准差，σ＝5.9；

(1)建立离散时间的粒子滤波模型为：

\{\begin{matrix} {x_{k}}^{(m)} = {x^{(m)}}_{k - 1} + {v^{(m)}}_{k} + {f_{k}}^{(m)} \\ {y_{k}}^{(m)} = h^{(m)} ({x_{k}}^{(m)}) + {ϵ^{(m)}}_{k} \end{matrix}

观测向量为

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T},

代表了测量得到的伪距与理论计算得到的伪距之间的差值；ε^(m) _k表示第k时刻的状态向量粒子的观测噪声；h^(m)(·)代表与之间的函数关系；

(2)计算未归一化的粒子权重的过程为：

其中，未归一化的粒子权重可由下式得到：

w_{k}^{(m)} = \frac{p (y_{k}^{(m)} | x_{k}^{(m)}) \hat{p} (x_{k}^{(m)} | x_{k - 1}^{(m)})}{q (x_{k}^{(m)} | x_{k - 1}^{(m)}, y_{k}^{(m)})} w_{k - 1}^{(m)}

w_{k}^{(m)} = p (y_{k}^{(m)} | x_{k}^{(m)}) w_{k - 1}^{(m)}

其中，经过重采样后，我们可以认为

w_{k - 1}^{(m)} = \frac{1}{M}

p (y_{k}^{(m)} | x_{k}^{(m)}) = {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}]

M＝10000不失一般性地假设第N颗卫星是欺骗卫星，从而得到了未归一化粒子权重的最终表达：

\begin{matrix} w_{k}^{(m)} = p (y_{k}^{(m)} | x_{k}^{(m)}) w_{k - 1}^{(m)} \\ = \frac{1}{M} {(\frac{1}{\sqrt{2 π σ^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}] \end{matrix}

归一化后为：

{\overset{&OverBar;}{w}}_{k}^{(m)} = \frac{w_{k}^{(m)}}{Σ_{m = 1}^{M} w_{k}^{(m)}};

\begin{matrix} θ = - 2 σ^{2} (\ln (M \max (w_{k}^{(m)})) + \frac{N}{2} \ln (2 π σ^{2})) \\ = \min {| | y_{k}^{(m)} | |}^{2} \end{matrix}

其中N为可见卫星的个数，取5；

{\hat{ρ}}_{f, k} = ρ_{f, k} + {Δρ}_{f, k}

\max (w_{k}^{(m)}) = \frac{1}{M} {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{\min (| | y_{k}^{(m)} | |)}{2 σ^{2}}]

这里的M为粒子数，取10000，根据定义监测变量θ为：

\begin{matrix} θ = - 2 σ^{2} (\ln (M \max (w_{k}^{(m)})) + \frac{N}{2} \ln (2 π σ^{2})) \\ = \min {| | y_{k}^{(m)} | |}^{2} \end{matrix}

当存在欺骗卫星时且可见卫星数为N＝5时，θ服从非中心卡方分布，其概率密度函数为：

g_{f} (θ) = \{\begin{matrix} \frac{1}{2 \sqrt{2 π σ^{2} θ}} \exp (- \frac{θ + λ}{2 σ^{2}}) \cdot (\exp (\frac{\sqrt{θλ}}{σ^{2}}) + \exp (- \frac{\sqrt{θλ}}{σ^{2}})), & θ &GreaterEqual; 0 \\ 0, & θ &GreaterEqual; 0 \end{matrix}

在N>5的条件下，g_f(θ)可以近似为：

g_{f} (θ) = \{\begin{matrix} \frac{1}{2 σ^{2}} {(\frac{θ}{λ})}^{\frac{N - 6}{4}} \exp (- \frac{θ + λ}{2 σ^{2}}) I_{\frac{N - 6}{2}} (\frac{\sqrt{θλ}}{σ^{2}}), & θ &GreaterEqual; 0 \\ 0, & θ &GreaterEqual; 0 \end{matrix}

I_T(·)为T阶修正的贝塞尔函数，

如果不存在欺骗卫星，则θ服从中心卡方分布，自由度为可见卫星数减4,即N-4：

g (θ) = \{\begin{matrix} \frac{1}{2 σ^{2} Γ (\frac{N}{2})} {(- \frac{θ}{{2 σ}^{2}})}^{\frac{N - 2}{2}} \exp (- \frac{θ}{2 σ^{2}}), & θ &GreaterEqual; 0 \\ 0, & θ < 0 \end{matrix}

λ = {(b_{k})}^{T} (I - G_{k} {(G_{k}^{T} G_{k})}^{- 1} G_{k}^{T}) b_{k}

b_k＝[0 0 ... b_s,k ... 0]^T

G_{k} = [\begin{matrix} 1_{1, k} & 1 \\ 1_{2, k} & 1 \\ . & . \\ . & . \\ . & . \\ 1_{N, k} & 1 \end{matrix}]

可见，存在欺骗卫星的情况下和不存在欺骗卫星的情况下，θ服从不同的分布，可以通过该分布来判断是否存在欺骗卫星；因此，可以设置一个门限γ，当θ>γ时认为存在欺骗卫星，反之则认为没有；这里的γ与可见卫星数N有关，N＝5时，γ取75，其他情况下γ取值可用试验方法得到；依据T阶修正的贝塞尔函数I_T(·)说明本方法的合理性如图3所示；

步骤七、判断θ与γ的大小，如果θ>γ，根据修正这里的γ取75，得到：

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix}

其中，b为误差下界取6，d为误差上界取8，为观测误差，P_f(·)为误差抑制函数；

P_{f} = {&Integral;}_{γ}^{+ \infty} g (θ) dθ,

相应的，发现概率为：

P_{D} = {&Integral;}_{γ}^{+ \infty} g_{f} (θ) dθ;

首先定义误差抑制函数P_f(·)：

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix}

经过这一处理后：

(1)若误差值的绝对值超过d(可以取d＝8，b＝6)，则说明对应的卫星是一颗假卫星，将它修正为0，即它就不会对后续计算产生任何影响了，也就是说不会影响到定位结果，所以假卫星被完全抑制，

(2)若的绝对值在b和d之间，则修正后减小为b，即

(3)若的绝对值小于b，则其大小则保持不变，即然后再重新计算未归一化粒子权重，得到归一化粒子权重这一部分的作用在于，抑制了欺骗卫星对定位结果的干扰；如果θ≤γ，则直接跳到步骤八；

步骤八、根据未归一化的粒子权重计算k时刻定位结果

考虑最小均方估计，则对于的估计为

{\overset{&OverBar;}{w}}_{k}^{(m)} = \frac{w_{k}^{(m)}}{Σ_{m = 1}^{M} w_{k}^{(m)}},

\begin{matrix} {\hat{x}}_{k}^{MMS} = Σ_{m = 1}^{M} {&Integral;}_{- \infty}^{+ \infty} y_{k}^{(m)} p (x_{k}^{(m)} | y_{k}^{(m)}) d x_{k}^{(m)} \\ = Σ_{m = 1}^{M} {\overset{&OverBar;}{w}}_{k}^{(m)} x_{k}^{(m)} \end{matrix};

其中，为归一化粒子权重；即作为第k个时刻的定位结果输出；

M_{eff} = \frac{1}{Σ_{m = 1}^{M} {(w_{k}^{(m)})}^{2}},

如果则转到步骤十一，如果跳到步骤十，进行重采样；

步骤十、如果采取M个在k时刻的状态向量粒子对集合重采样，采到的概率为归一化粒子权重重新采到的归一化粒子权重直至使得为止的具体过程为：

对集合中的元素进行均匀采样，依次抽取M次，得到k时刻M个新的状态向量粒子，用k时刻M个新的状态向量粒子来取代之前的k时刻的M个状态向量粒子其中M个状态向量粒子抽取的次数为M次，如果则说明有效粒子数M_eff不足，需要进行重采样，以解决粒子退化问题；所谓粒子退化问题是指，经过若干次递推后，只有少数粒子的权重很大，而大部分粒子的权重非常接近于零的现象；由于权重接近于零的粒子对状态的估计的贡献很小，所以当粒子退化严重时，有很多计算和存储资源都被浪费了；

步骤十一、如果利用步骤二中的概率为的在k时刻状态向量粒子将k时刻加1后得到用户在k+1时刻的用户状态向量粒子跳到步骤二，开始下一个时刻的定位过程；

仿真结果表明，该发明可以有效侦测和抑制欺骗卫星的干扰,首先给出侦测效果的仿真结果,仿真参数设置如下表：

变量	取值
		σ	5.9m
M	10000
		计算频率	1次/秒
仿真时间	250s

实施例二

x_{1}^{(m)} = x_{0}^{(m)} + v_{1}^{(m)} + f_{1}^{(m)},

其中，

v_{1} = {[{\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, {\overset{\cdot}{δ} t}_{U, 1}]}^{T}, {\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, \overset{\cdot}{δ} t_{U, 1}

x_{k}^{(m)} = x_{k - 1}^{(m)} + v_{k}^{(m)} + f_{k}^{(m)}; k = 2,3 . . . . .

(1)通过用户在第k时刻的状态向量粒子计算得到观测伪距

步骤四、根据步骤三得到的观测伪距更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T}

步骤五、根据更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T}

计算未归一化的粒子权重

w_{k}^{(m)} = \frac{1}{M} {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}]

(1)建立离散时间的粒子滤波模型为：

\{\begin{matrix} {x_{k}}^{(m)} = {x^{(m)}}_{k - 1} + {v^{(m)}}_{k} + {f_{k}}^{(m)} \\ {y_{k}}^{(m)} = h^{(m)} ({x_{k}}^{(m)}) + {ϵ^{(m)}}_{k} \end{matrix}

(2)计算未归一化的粒子权重的过程为：

其中，未归一化的粒子权重可由下式得到：

w_{k}^{(m)} = \frac{p (y_{k}^{(m)} | x_{k}^{(m)}) \hat{p} (x_{k}^{(m)} | x_{k - 1}^{(m)})}{q (x_{k}^{(m)} | x_{k - 1}^{(m)}, y_{k}^{(m)})} w_{k - 1}^{(m)}

w_{k}^{(m)} = p (y_{k}^{(m)} | x_{k}^{(m)}) w_{k - 1}^{(m)}

其中，经过重采样后，我们可以认为

w_{k - 1}^{(m)} = \frac{1}{M}

p (y_{k}^{(m)} | x_{k}^{(m)}) = {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}]

\begin{matrix} w_{k}^{(m)} = p (y_{k}^{(m)} | x_{k}^{(m)}) w_{k - 1}^{(m)} \\ = \frac{1}{M} {(\frac{1}{\sqrt{2 π σ^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}] \end{matrix}

其中，M为在第k时刻的状态向量粒子数，取10000，N为可见卫星的个数，取8，σ为的标准差，σ＝5.9；

归一化后为：

{\overset{&OverBar;}{w}}_{k}^{(m)} = \frac{w_{k}^{(m)}}{Σ_{m = 1}^{M} w_{k}^{(m)}};

\begin{matrix} θ = - 2 σ^{2} (\ln (M \max (w_{k}^{(m)})) + \frac{N}{2} \ln (2 π σ^{2})) \\ = \min {| | y_{k}^{(m)} | |}^{2} \end{matrix}

其中N为可见卫星的个数，取5；

{\hat{ρ}}_{f, k} = ρ_{f, k} + {Δρ}_{f, k}

\max (w_{k}^{(m)}) = \frac{1}{M} {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{\min (| | y_{k}^{(m)} | |)}{2 σ^{2}}]

这里的M为粒子数，取10000，根据定义监测变量θ为：

\begin{matrix} θ = - 2 σ^{2} (\ln (M \max (w_{k}^{(m)})) + \frac{N}{2} \ln (2 π σ^{2})) \\ = \min {| | y_{k}^{(m)} | |}^{2} \end{matrix}

g_{f} (θ) = \{\begin{matrix} \frac{1}{2 \sqrt{2 π σ^{2} θ}} \exp (- \frac{θ + λ}{2 σ^{2}}) \cdot (\exp (\frac{\sqrt{θλ}}{σ^{2}}) + \exp (- \frac{\sqrt{θλ}}{σ^{2}})), & θ &GreaterEqual; 0 \\ 0, & θ &GreaterEqual; 0 \end{matrix}

在N>5的条件下，g_f(θ)可以近似为：

g_{f} (θ) = \{\begin{matrix} \frac{1}{2 σ^{2}} {(\frac{θ}{λ})}^{\frac{N - 6}{4}} \exp (- \frac{θ + λ}{2 σ^{2}}) I_{\frac{N - 6}{2}} (\frac{\sqrt{θλ}}{σ^{2}}), & θ &GreaterEqual; 0 \\ 0, & θ &GreaterEqual; 0 \end{matrix}

I_T(·)为T阶修正的贝塞尔函数，

g (θ) = \{\begin{matrix} \frac{1}{2 σ^{2} Γ (\frac{N}{2})} {(- \frac{θ}{{2 σ}^{2}})}^{\frac{N - 2}{2}} \exp (- \frac{θ}{2 σ^{2}}), & θ &GreaterEqual; 0 \\ 0, & θ < 0 \end{matrix}

λ = {(b_{k})}^{T} (I - G_{k} {(G_{k}^{T} G_{k})}^{- 1} G_{k}^{T}) b_{k}

b_k＝[0 0 ... b_s,k ... 0]^T

G_{k} = [\begin{matrix} 1_{1, k} & 1 \\ 1_{2, k} & 1 \\ . & . \\ . & . \\ . & . \\ 1_{N, k} & 1 \end{matrix}]

可见，存在欺骗卫星的情况下和不存在欺骗卫星的情况下，θ服从不同的分布，可以通过该分布来判断是否存在欺骗卫星；因此，可以设置一个门限γ，当θ>γ时认为存在欺骗卫星，反之则认为没有；这里的γ与可见卫星数N有关，N＝8时，γ取250，其他情况下γ取值可用试验方法得到；依据T阶修正的贝塞尔函数I_T(·)说明本方法的合理性如图4所示；

步骤七、判断θ与γ的大小，如果θ>γ，根据修正这里的γ取250，得到：

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix}

P_{f} = {&Integral;}_{γ}^{+ \infty} g (θ) dθ,

相应的，发现概率为：

P_{D} = {&Integral;}_{γ}^{+ \infty} g_{f} (θ) dθ;

首先定义误差抑制函数P_f(·)：

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix}

经过这一处理后：

(2)若的绝对值在b和d之间，则修正后减小为b，即

步骤八、根据未归一化的粒子权重计算k时刻定位结果

考虑最小均方估计，则对于的估计为

{\overset{&OverBar;}{w}}_{k}^{(m)} = \frac{w_{k}^{(m)}}{Σ_{m = 1}^{M} w_{k}^{(m)}},

\begin{matrix} {\hat{x}}_{k}^{MMS} = Σ_{m = 1}^{M} {&Integral;}_{- \infty}^{+ \infty} y_{k}^{(m)} p (x_{k}^{(m)} | y_{k}^{(m)}) d x_{k}^{(m)} \\ = Σ_{m = 1}^{M} {\overset{&OverBar;}{w}}_{k}^{(m)} x_{k}^{(m)} \end{matrix};

M_{eff} = \frac{1}{Σ_{m = 1}^{M} {(w_{k}^{(m)})}^{2}},

如果则转到步骤十一，如果跳到步骤十，进行重采样；

对集合中的元素进行均匀采样，依次抽取M次，得到k时刻M个新的状态向量粒子，用k时刻M个新的状态向量粒子来取代之前的k时刻的M个状态向量粒子其中M个状态向量粒子抽取的次数为M次；如果则说明有效粒子数M_eff不足，需要进行重采样，以解决粒子退化问题；所谓粒子退化问题是指，经过若干次递推后，只有少数粒子的权重很大，而大部分粒子的权重非常接近于零的现象；由于权重接近于零的粒子对状态的估计的贡献很小，所以当粒子退化严重时，有很多计算和存储资源都被浪费了；

变量	取值
		σ	5.9m
M	10000
		计算频率	1次/秒
仿真时间	250s

Claims

1.一种基于粒子滤波的稳健GNSS抗欺骗方法，其特征在于：一种基于粒子滤波的稳健GNSS抗欺骗方法具体是按照以下步骤进行的：

x_{1}^{(m)} = x_{0}^{(m)} + v_{1}^{(m)} + f_{1}^{(m)},

其中，

v_{1} = {[{\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, {\overset{\cdot}{δ} t}_{U, 1}]}^{T}, {\overset{\cdot}{x}}_{U, 1}, {\overset{\cdot}{y}}_{{U,}_{1}}, {\overset{\cdot}{z}}_{U, 1}, \overset{\cdot}{δ} t_{U, 1}

x_{k}^{(m)} = x_{k - 1}^{(m)} + v_{k}^{(m)} + f_{k}^{(m)}; k = 2,3 . . . . .

其中，状态向量粒子(x_U,k,y_U,k,z_U,k)代表用户在k个时刻状态粒子的三维位置；δt_U,k用户在k个时刻状态与第n颗卫星的时钟偏差；模型输入为：

v_{k}^{(m)} = {[{\overset{\cdot}{x}}_{U, k}, {\overset{\cdot}{y}}_{U, k}, {\overset{\cdot}{z}}_{U, k}, {\overset{\cdot}{δ} t}_{U, k}]}^{T}, {\overset{\cdot}{x}}_{U, k}, {\overset{\cdot}{y}}_{U, k}, {\overset{\cdot}{z}}_{U, k}, \overset{\cdot}{δ} t_{U, k}

步骤四、根据步骤三得到的观测伪距更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T}

步骤五、根据更新观测向量

y_{k}^{(m)} = {[e_{1, k}^{(m)}, . . ., e_{n, k}^{(m)}, . . ., e_{N, k}^{(m)}]}^{T}

计算未归一化的粒子权重

w_{k}^{(m)} = \frac{1}{M} {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}]

\begin{matrix} θ = - 2 σ^{2} (\ln (M \max (w_{k}^{(m)})) + \frac{N}{2} \ln (2 π σ^{2})) \\ = \min {| | y_{k}^{(m)} | |}^{2} \end{matrix};

步骤七、判断θ与γ的大小，如果θ>γ，根据修正得到：

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix}

然后根据修正后的再重新计算未归一化的粒子权重

如果θ≤γ，则直接跳到步骤八；

步骤八、根据未归一化的粒子权重计算k时刻定位结果

{\overset{&OverBar;}{w}}_{k}^{(m)} = \frac{w_{k}^{(m)}}{Σ_{m = 1}^{M} w_{k}^{(m)}},

{\hat{x}}_{k}^{MMS} = Σ_{m = 1}^{M} {\overset{&OverBar;}{w}}_{k}^{(m)} x_{k}^{(m)};

其中，为归一化粒子权重；

M_{eff} = \frac{1}{Σ_{m = 1}^{M} {(w_{k}^{(m)})}^{2}},

如果则转到步骤十一，如果跳到步骤十，进行重采样；

2.根据权利要求1所述一种基于粒子滤波的稳健GNSS抗欺骗方法，其特征在于：步骤三中根据步骤二计算的计算出用户与第n颗卫星之间的观测伪距由下式给出：

(1)通过用户在第k时刻的状态向量粒子计算得到观测伪距

其中，(x_n,k,y_n,k,z_n,k)(n＝1,2,...,N)代表在第n颗卫星在k个时刻的位置坐标；(x_U,k,y_U,k,z_U,k)为在第k时刻的状态粒子三维位置，δt_U,k为相对于第n颗卫星的时钟偏差；

(2)真实伪距ρ_n,k由接收机的伪距测量电路直接测量得到，假设可见卫星的数目为N，在第k个时刻，测得用户和第n颗卫星之间的观测伪距与真实伪距ρ_n,k的关系为：

ε_n,k表示在k时刻接收机测量与第n颗卫星的伪距时产生的测量误差。

3.根据权利要求1所述一种基于粒子滤波的稳健GNSS抗欺骗方法，其特征在于：步骤五中根据更新观测向量计算未归一化的粒子权重的过程为：

(1)建立离散时间的粒子滤波模型为：

\{\begin{matrix} {x_{k}}^{(m)} = {x^{(m)}}_{k - 1} + {v^{(m)}}_{k} + {f_{k}}^{(m)} \\ {y_{k}}^{(m)} = h^{(m)} ({x_{k}}^{(m)}) + {ϵ^{(m)}}_{k} \end{matrix}

(2)计算未归一化的粒子权重的过程为：

其中，未归一化的粒子权重可由下式得到：

w_{k}^{(m)} = \frac{p (y_{k}^{(m)} | x_{k}^{(m)}) \hat{p} (x_{k}^{(m)} | x_{k - 1}^{(m)})}{q (x_{k}^{(m)} | x_{k - 1}^{(m)}, y_{k}^{(m)})} w_{k - 1}^{(m)}

上式中，q(·)为重要密度函数，每一个状态向量粒子都来自对q(·)的完备采样，为已知时的概率密度函数，为已知时的概率密度函数，得到了未归一化粒子权重的表达式：

w_{k}^{(m)} = p (y_{k}^{(m)} | x_{k}^{(m)}) w_{k - 1}^{(m)}

其中，经过重采样后：

w_{k - 1}^{(m)} = \frac{1}{M}

p (y_{k}^{(m)} | x_{k}^{(m)}) = {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}]

从而得到了未归一化粒子权重的最终表达：

\begin{matrix} w_{k}^{(m)} = p (y_{k}^{(m)} | x_{k}^{(m)}) w_{k - 1}^{(m)} \\ = \frac{1}{M} {(\frac{1}{\sqrt{2 π σ^{2}}})}^{N} \exp [- \frac{| | y_{k}^{(m)} | |}{2 σ^{2}}] \end{matrix}

归一化后为：

{\overset{&OverBar;}{w}}_{k}^{(m)} = \frac{w_{k}^{(m)}}{Σ_{m = 1}^{M} w_{k}^{(m)}} .

4.根据权利要求1所述一种基于粒子滤波的稳健GNSS抗欺骗方法，其特征在于：步骤六中计算观测向量的模值的最小值θ：

通过监测来侦测欺骗卫星，具体如下：

\max (w_{k}^{(m)}) = \frac{1}{M} {(\frac{1}{\sqrt{2 {πσ}^{2}}})}^{N} \exp [- \frac{\min (| | y_{k}^{(m)} | |)}{2 σ^{2}}]

根据定义监测变量θ为：

\begin{matrix} θ = - 2 σ^{2} (\ln (M \max (w_{k}^{(m)})) + \frac{N}{2} \ln (2 π σ^{2})) \\ = \min {| | y_{k}^{(m)} | |}^{2} \end{matrix}

g_{f} (θ) = \{\begin{matrix} \frac{1}{2 \sqrt{2 π σ^{2} θ}} \exp (- \frac{θ + λ}{2 σ^{2}}) \cdot (\exp (\frac{\sqrt{θλ}}{σ^{2}}) + \exp (- \frac{\sqrt{θλ}}{σ^{2}})), & θ &GreaterEqual; 0 \\ 0, & θ &GreaterEqual; 0 \end{matrix}

在N>5的条件下，g_f(θ)可以近似为：

g_{f} (θ) = \{\begin{matrix} \frac{1}{2 σ^{2}} {(\frac{θ}{λ})}^{\frac{N - 6}{4}} \exp (- \frac{θ + λ}{2 σ^{2}}) I_{\frac{N - 6}{2}} (\frac{\sqrt{θλ}}{σ^{2}}), & θ &GreaterEqual; 0 \\ 0, & θ &GreaterEqual; 0 \end{matrix}

I_T(·)为T阶修正的贝塞尔函数；

g (θ) = \{\begin{matrix} \frac{1}{2 σ^{2} Γ (\frac{N}{2})} {(- \frac{θ}{{2 σ}^{2}})}^{\frac{N - 2}{2}} \exp (- \frac{θ}{2 σ^{2}}), & θ &GreaterEqual; 0 \\ 0, & θ < 0 \end{matrix}

其中，Γ(·)是gamma函数，λ为卡方分布的非中心度参数，其取值与可见卫星个数N有关，也与伪距有关，由下面的方法得出：

λ = {(b_{k})}^{T} (I - G_{k} {(G_{k}^{T} G_{k})}^{- 1} G_{k}^{T}) b_{k}

b_k＝[0 0 ... b_s,k ... 0]^T

G_{k} = [\begin{matrix} 1_{1, k} & 1 \\ 1_{2, k} & 1 \\ . & . \\ . & . \\ . & . \\ 1_{N, k} & 1 \end{matrix}]

在第k个时刻，这里1_n,k是当前接收机指向第n颗卫星的单位方向矢量，I表示单位矩阵，b_s,k为第n颗卫星的附加伪距，b_k为对应的附加伪距矩阵，G_k表示接收机与所有可见卫星的单位方向矩阵，T代表矩阵的转置，是G_k的转置。

5.根据权利要求1所述一种基于粒子滤波的稳健GNSS抗欺骗方法，其特征在于：步骤七中判断θ与γ的大小，如果θ>γ，根据修正得到

其中

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix},

P_f(·)为误差抑制函数，为观测误差，b为误差下界，d为误差上界；

然后根据修正后的再重新计算未归一化粒子权重

虚警概率P_f与门限γ的关系由下面的式子给出：

P_{f} = {&Integral;}_{γ}^{+ \infty} g (θ) dθ,

相应的，发现概率为：

P_{D} = {&Integral;}_{γ}^{+ \infty} g_{f} (θ) dθ;

在侦测到欺骗卫星后，可以通过基于稳健估计的改进方法来对其进行抑制具体处理方法如下：

首先定义误差抑制函数P_f(·)：

P_{f} (e_{n, k}^{(m)}) = \{\begin{matrix} 1, & b > | e_{n, k}^{(m)} | \\ \frac{b}{| e_{n, k}^{(m)} |}, & b > | e_{n, k}^{(m)} | > d \\ 0, & | e_{n, k}^{(m)} | > d \end{matrix}

其中，b为误差下界，d为误差上界，在计算得到未归一化粒子权重后，对观测向量中的误差值进行修正，得到方法如下：

经过这一处理后：

(1)若误差值的绝对值超过d则说明对应的卫星是一颗假卫星，将它修正为0，即

(2)若的绝对值在b和d之间，则修正后减小为b，即

(3)若的绝对值小于b，则其大小则保持不变，即然后再重新计算未归一化粒子权重，得到归一化粒子权重

6.根据权利要求1所述一种基于粒子滤波的稳健GNSS抗欺骗方法，其特征在于：步骤八中根据未归一化的粒子权重计算k时刻定位结果推导过程为：

则对于的估计为

\begin{matrix} {\hat{x}}_{k}^{MMS} = Σ_{m = 1}^{M} {&Integral;}_{- \infty}^{+ \infty} y_{k}^{(m)} p (x_{k}^{(m)} | y_{k}^{(m)}) d x_{k}^{(m)} \\ = Σ_{m = 1}^{M} {\overset{&OverBar;}{w}}_{k}^{(m)} x_{k}^{(m)} \end{matrix},

即作为第k个时刻的定位结果输出。

7.根据权利要求1所述一种基于粒子滤波的稳健GNSS抗欺骗方法，其特征在于：步骤十中如果采取M个在k时刻的状态向量粒子对集合重采样，采到的概率为归一化粒子权重具体过程为：

对集合中的元素进行均匀采样，依次抽取M次，得到k时刻M个新的状态向量粒子，用k时刻M个新的状态向量粒子来取代之前的k时刻的M个状态向量粒子