法律状态公告日
法律状态信息
法律状态
2016-05-18
授权
授权
2014-08-13
实质审查的生效 IPC(主分类):G01S19/21 申请日:20140425
实质审查的生效
2014-07-16
公开
公开
技术领域
本发明涉及一种基于粒子滤波的稳健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时刻用户状态向量粒子
其中,
步骤二、当k≠1时,根据k-1时刻的用户状态向量粒子更新m个在第k时刻状 态向量粒子
其中,状态向量粒子(xU,k,yU,k,zU,k)代表用户在k个时刻状态向量粒子的三维 位置;δtU,k用户在k个时刻状态与第n颗卫星的时钟偏差;模型输入为:
步骤三、根据步骤二计算的计算出用户与第n颗卫星之间的观测伪距由 下式给出:
是接收机和卫星之间的观测伪距,(xn,k,yn,k,zn,k)(n=1,2,...,N)代表第 n颗卫星在第k个时刻的位置坐标,c表示光速;
步骤四、根据步骤三得到的观测伪距更新观测向量
其中,用户与第n颗卫星的伪距观测误差为ρn,k表示在第k个时 刻,第n颗卫星与用户之间的真实伪距,N为可见卫星的个数;
步骤五、根据更新观测向量
计算未归一化的粒子权重
其中,M为在第k时刻的状态向量粒子数,N为可见卫星的个数,σ为的标准差;
步骤六、根据未归一化的粒子权重计算观测向量的模值的最小值θ:
步骤七、判断θ与γ的大小,如果θ>γ,根据修正得到:
其中,b为误差下界,d为误差上界,为观测误差,Pf(·)为误差抑制函数;
然后根据修正后的再重新计算未归一化的粒子权重如果θ≤γ,则直接跳到步骤八;
步骤八、根据未归一化的粒子权重计算k时刻定位结果
其中,为归一化粒子权重;
步骤九、将步骤二中共得到M个粒子在k时刻的状态向量粒子根据步骤五中的 未归一化的粒子权重计算得到剩余的有效向量粒子的数目Meff:
如果则转到步骤十一,如果跳到步骤十,进行重采样;
步骤十、如果采取M个在k时刻的状态向量粒子对集合重采样,采到的概率为归一化粒子权重
步骤十一、如果利用步骤二中的概率为的在k时刻状态向量粒子 将k时刻加1后得到用户在k+1时刻的用户状态向量粒子跳到步骤二,开 始下一个时刻的定位过程;即完成了一种基于粒子滤波的稳健GNSS抗欺骗方法。
发明效果
本发明提出了一种基于粒子滤波的抗欺骗方法,能够在不增加任何硬件资源、不改变 信号体制的条件下实现欺骗信号的侦查和消除。从而确保了定位精度的准确,保护了国家 和人民的利益。该发明可以有效侦测和抑制欺骗卫星的干扰,图2可以看到,可见卫星越 多、欺骗卫星给出的伪距误差越大,则检测效果越好。在误差达到40m时,检测概率可 以达到70%以上。为了证明该算法对欺骗卫星干扰的抑制,下面对传统的最小二乘法、 传统的粒子滤波算法和本文的算法做了比较。图5、图6可以看到,仿真时间为250秒, 前120秒没有干扰,之后加入一颗欺骗卫星,仿真时间为250秒0秒到120秒没有干扰, 三种方法的定位误差差不多,在第120秒加入干扰后,传统的最小二乘法和粒子滤波方法 都产生了很大的定位误差,但本文的算法却还是保持了很好的定位精度。转发式的干扰, 其主要机理是提供假的伪距,从而欺骗接收机。而本发明的本质,就是检测伪距的残差, 通过伪距的异常来判断是否存在欺骗卫星。可见,本发明可以检测一切提供假伪距的欺骗 方法,转发式干扰自然也包括在内。本发明不需要新的信号架构,即不需要改变信号体制, 可以直接应用于已有的系统。这一点从后面的分析过程中可以看出,该方法使用的参数都 是现有系统可以提供的
本发明直接应用于已有系统,不需要任何额外的硬件,只需要修改解算过程,即对接 收机的软件做调整即可。而现在多数其他的抗欺骗式干扰的方法是需要额外的硬件的,所 以本发明相对来说成本低。
附图说明
图1是具体实施方式一提出的一种基于粒子滤波的稳健GNSS抗欺骗方法流程图;
图2是具体实施方式一提出的发现欺骗卫星的概率PD与卫星数N的关系图,其中横 坐标△Pf为欺骗卫星提供的伪距和卫星真实伪距的偏差,纵坐标为发现欺骗卫星的概率 PD;为自由度为N=5时的发现欺骗卫星的概率PD,为自由度为N=8时的发现 欺骗卫星的概率PD;
图3是具体实施方式四提出的自由度N=5时g(θ)与gf(θ)的关系图;其中,横坐标为 监测变量θ,纵坐标PDF of θ为监测变量θ的概率密度函数值;为不存在欺骗卫星, 监测变量θ的概率密度函数g(θ)值变化曲线;为存在欺骗卫星,则监测变量θ的概 率密度函数gf(θ)值的变化曲线;
图4是具体实施方式四提出的N=8时g(θ)与gf(θ)的关系图;其中,横坐标为监测变 量θ,纵坐标PDF of θ为监测变量θ的概率密度函数值;为不存在欺骗卫星,监测 变量θ的概率密度函数g(θ)值变化曲线;为存在欺骗卫星,则监测变量θ的概率密 度函数gf(θ)值的变化曲线;
图5是具体实施方式一提出的N=5时不同算法的定位误差与定位时间曲线图;其中 为传统的最小二乘法,为传统的抗干扰方法,为本发明提出的抗干扰 算法;
图6是具体实施方式一提出的N=8时不同算法的定位误差与定位时间曲线图;其中 为传统的最小二乘法,为传统的抗干扰方法,为本发明提出的抗干扰 算法。
具体实施方式
具体实施方式一:本实施方式的一种基于粒子滤波的稳健GNSS抗欺骗方法,具体 是按照以下步骤制备的:
步骤一、完成两个变量的初始化,首先将确定用户的初始状态,在用户初始状态的领 域内均匀分布m个初始状态向量粒子均匀分布在x0领域内;同 时设置迭代次数的计数值为k=1;根据用户初始状态向量粒子可以得到k=1时刻用 户状态向量粒子
其中,
步骤二、当k≠1时,根据k-1时刻的用户状态向量粒子更新m个在第k时刻状 态向量粒子
是粒子滤波的动态处理模型,以适应GNSS接收机的运动;其中,状态向量粒子 也就是用户在k时隙的状态,也是最终的定位结果;所有坐 标都是在ECEF(Earth-Centered and Earth-Fixed,地心固定坐标系)坐标系下给出, (xU,k,yU,k,zU,k)代表用户在k个时刻状态粒子的三维位置;δtU,k用户在k个时刻状态 与第n颗卫星的时钟偏差;模型输入为:分别表示xU,k,yU,k,zU,k,δtU,k在单位时隙中的变化量;表示用 户在第k时刻的状态向量粒子的过程噪声,服从均值为零的高斯分布;
步骤三、根据步骤二计算的计算出用户与第n颗卫星之间的观测伪距由 下式给出:
是接收机和卫星之间的观测伪距,(xn,k,yn,k,zn,k)(n=1,2,...,N)代表第 n颗卫星在第k个时刻的位置坐标,c表示光速;
步骤四、根据步骤三得到的观测伪距更新观测向量
其中,用户与第n颗卫星的伪距观测误差为ρn,k表示在第k个时 刻,第n颗卫星与用户之间的真实伪距,N为可见卫星的个数;
步骤五、根据更新观测向量
计算未归一化的粒子权重
其中,M为在第k时刻的状态向量粒子数,N为可见卫星的个数,σ为的标准差;
步骤六、根据未归一化的粒子权重计算观测向量的模值的最小值θ:
步骤七、判断θ与γ的大小,如果θ>γ,根据修正得 到:
其中
其中,b为误差下界,d为误差上界,为观测误差,Pf(·)为误差抑制函数;
然后根据修正后的再重新计算未归一化的粒子权重使得θ≤γ;
如果θ≤γ,则直接跳到步骤八;
步骤八、根据未归一化的粒子权重计算k时刻定位结果
其中,为归一化粒子权重;
步骤九、将步骤二中共得到M个粒子在k时刻的状态向量粒子由于存在粒子退 化现象,部分状态粒子将失效,剩余的有效状态粒子的数目Meff,根据步骤五中的未归 一化的粒子权重计算得到剩余的有效向量粒子的数目Meff:
如果则转到步骤十一,如果跳到步骤十,进行重采样;
步骤十、如果采取M个在k时刻的状态向量粒子对集合 重采样,采到的概率为归一化粒子权重重新采到的归一化粒子权重直至使得 为止;
步骤十一、如果利用步骤二中的概率为的在k时刻状态向量粒子 将k时刻加1后得到用户在k+1时刻的用户状态向量粒子跳到步骤二,开 始下一个时刻的定位过程;如图1即完成了一种基于粒子滤波的稳健GNSS抗欺骗方法。
本实施方式效果:
本实施方式提出了一种基于粒子滤波的抗欺骗方法,能够在不增加任何硬件资源、不 改变信号体制的条件下实现欺骗信号的侦查和消除。从而确保了定位精度的准确,保护了 国家和人民的利益。本实施方式可以有效侦测和抑制欺骗卫星的干扰,图2可以看到,可 见卫星越多、欺骗卫星给出的伪距误差越大,则检测效果越好。在误差达到40m时,检 测概率可以达到70%以上。为了证明该算法对欺骗卫星干扰的抑制,下面对传统的最小 二乘法、传统的粒子滤波算法和本文的算法做了比较。图5、图6可以看到,仿真时间为 250秒,前120秒没有干扰,之后加入一颗欺骗卫星,仿真时间为250秒0秒到120秒没 有干扰,三种方法的定位误差差不多,在第120秒加入干扰后,传统的最小二乘法和粒子 滤波方法都产生了很大的定位误差,但本文的算法却还是保持了很好的定位精度。转发式 的干扰,其主要机理是提供假的伪距,从而欺骗接收机。而本实施方式的本质,就是检测 伪距的残差,通过伪距的异常来判断是否存在欺骗卫星。可见,本实施方式可以检测一切 提供假伪距的欺骗方法,转发式干扰自然也包括在内。本实施方式不需要新的信号架构, 即不需要改变信号体制,可以直接应用于已有的系统。这一点从后面的分析过程中可以看 出,该方法使用的参数都是现有系统可以提供的
本实施方式直接应用于已有系统,不需要任何额外的硬件,只需要修改解算过程,即 对接收机的软件做调整即可。而现在多数其他的抗欺骗式干扰的方法是需要额外的硬件 的,所以本发明相对来说成本低。
具体实施方式二:本实施方式与具体实施方式一不同的是:步骤三中根据步骤二计算 的计算出用户与第n颗卫星之间的观测伪距由下式给出:
(1)通过用户在第k时刻的状态向量粒子计算得到观测伪距
其中,(xn,k,yn,k,zn,k)(n=1,2,...,N)代表在第n颗卫星在k个时刻的 位置坐标;(xU,k,yU,k,zU,k,δtU,k)代表用户在第k时刻的状态向量粒子的状态, 包括用户在第k时刻的状态向量粒子三维位置(xU,k,yU,k,zU,k)和相对于第n颗卫星 的时钟偏差δtU,k;
(2)真实伪距ρn,k由接收机的伪距测量电路直接测量得到,假设可见卫星的数目为N, 信号捕获和跟踪的影响被忽略;在第k个时刻,测得用户和第n颗卫星之间的观测伪距与真实伪距ρn,k的关系为:
εn,k表示在k时刻接收机测量与第n颗卫星的伪距时产生的测量误差。其它步骤及参 数与具体实施方式一相同。
具体实施方式三:本实施方式与具体实施方式一或二不同的是:步骤五中根据更新观 测向量计算未归一化的粒子权重的过程为:
(1)建立离散时间的粒子滤波模型为:
观测向量为代表了测量得到的伪距与理论计 算得到的伪距之间的差值;ε(m)k表示第k时刻的状态向量粒子的观测噪声;h(m)(·) 代表与yk(m)之间的函数关系;
(2)计算未归一化的粒子权重的过程为:
其中,未归一化的粒子权重可由下式得到:
上式中,q(·)为重要密度函数(这个是为了推导方便,人为定义的一个函数,没有明 确的意义),每一个状态向量粒子都来自对q(·)的完备采样,为已知 时的概率密度函数,为已知时的概率密度函数;为了 得到方差最小的未归一化粒子权重,应该选取的后验概率密度 即已知和时的概率密度,但在大多数实际情况中是 无法实现的,可以采用近似表示,得到了未归一化粒子权重的表达式:
其中,经过重采样后,我们可以认为
不失一般性地假设第N颗卫星是欺骗卫星,从而得到了未归一化粒子权重的最 终表达:
归一化后为:
其它步骤及参数与具体实施方式一或二相同。
具体实施方式四:本实施方式与具体实施方式一至三之一不同的是:步骤六中计算观 测向量的模值的最小值θ:
假设有多于四颗可见的真实卫星,还有一颗欺骗卫星;这颗欺骗卫星会给接收机提供 一个自身坐标(xf,k,yf,k,zf,k)和假的伪距从而提出了对欺骗的侦测和抑制方法, 由下式给出:
ρf,k是接收机和卫星之间的真实伪距,是欺骗卫星提供的假伪距,△ρf,k是二者的偏 差;一般而言,定位误差与|△ρf,k|正相关;(本公式是理论上的推导,说明本方法的合理 性)
前面已经提到,就是测量伪距与理论计算得到的伪距之间的差值;在没有欺骗卫 星的条件下,该向量中的每一个数都不会很大;如果存在欺骗卫星,对应的伪距误差就会 很大,从而导致的二范数很大,也就导致很大;所以,通过监测来侦测欺 骗卫星,具体如下:
根据定义监测变量θ为:
当存在欺骗卫星时且自由度N=5时,θ服从非中心卡方分布,其概率密度函数为:
在N>5的条件下,gf(θ)可以近似为:
IT(·)为T阶修正的贝塞尔函数,
如果不存在欺骗卫星,则θ服从中心卡方分布,自由度为N-4:
其中,Γ(·)是gamma函数,λ为卡方分布的非中心度参数,其取值与可见卫星个数 N有关,也与伪距有关;由下面的方法得出:
bk=[0 0 ... bs,k ... 0]T
在第k个时刻,这里1n,k是当前接收机指向第n颗卫星的单位方向矢量,I表示单位矩 阵,bs,k为第n颗卫星的附加伪距,bk为对应的附加伪距矩阵,Gk表示接收机与所有可 见卫星的单位方向矩阵,T代表矩阵的转置,是Gk的转置;
可见,存在欺骗卫星的情况下和不存在欺骗卫星的情况下,θ服从不同的分布,可以 通过该分布来判断是否存在欺骗卫星;因此,可以设置一个门限γ,当θ>γ时认为存在 欺骗卫星,反之则认为没有;依据T阶修正的贝塞尔函数IT(·)说明本方法的合理性如图 3、4所示。其它步骤及参数与具体实施方式一至三之一相同。
具体实施方式五:本实施方式与具体实施方式一至四之一不同的是:步骤七中判断θ 与γ的大小,如果θ>γ,根据修正得到
其中
然后根据修正后的再重新计算未归一化粒子权重
当不存在欺骗卫星时,该方法有一定的概率发生误判,这一误判的概率被称为虚警概 率Pf,与门限γ有关;虚警概率Pf与门限γ的关系由下面的式子给出:
相应的,发现概率为:
其中,g(θ):不存在欺骗卫星的情况下,监测变量θ的概率密度函数;PD:发现欺 骗卫星的概率;gf(θ):存在欺骗卫星的情况下,监测变量θ的概率密度函数;
在侦测到欺骗卫星后,可以通过基于稳健估计的改进方法来对其进行抑制;稳健估计 的目的在于,在一些普遍假设的条件下,给出最贴近真实值的估计;传统的估计方法,比 如最小二乘法,认为每一个观测数据都服从先验条件,所以,某一个错误的观测值会带来 很大的误差;而稳健估计则不同,它会去除某些异常观测值,从而消除其影响;最大似然 类型的估计(M-estimate),高阶统计的线性组合(L-estimate),等,是基本的稳健估计方法; 其中,最大似然比估计总是被用在参数估计上;
本发明基于传统的稳健估计方法,提出了改进,具体处理方法如下:
首先定义误差抑制函数Pf(·):
其中,b为误差下界,d为误差上界,e为观测误差,在计算得到未归一化粒子权重 后,对观测误差进行修正,得到方法如下:
经过这一处理后:
(1)若误差值的绝对值超过d,则说明对应的卫星是一颗假卫星,将它修正为 0,即它就不会对后续计算产生任何影响了,也就是说 不会影响到定位结果,所以假卫星被完全抑制,
(2)若的绝对值在b和d之间,则修正后减小为b,即
(3)若的绝对值小于b,则其大小则保持不变,即 然后再重新计算未归一化粒子权重,得到归一化粒 子权重这一部分的作用在于,抑制了欺骗卫星对定位结果的干扰。其它步骤及参 数与具体实施方式一至四之一相同。
具体实施方式六:本实施方式与具体实施方式一至五之一不同的是:步骤八中根据未 归一化的粒子权重计算k时刻定位结果推导过程为:
考虑最小均方估计,则对于的估计为
即作为第k个时刻的定位结果输出。其它步骤及参数与具体实施方式一至五之 一相同。
具体实施方式七:本实施方式与具体实施方式一至六之一不同的是:步骤十中如果 采取M个在k时刻的用户状态向量粒子具体过程为:
对集合中的元素进行均匀采样,依次抽取M次,得到k时刻M个新的状态向 量粒子,用k时刻M个新的状态向量粒子来取代之前的k时刻的M个状态向量粒子其中M个状态向量粒子抽取的次数为M次,如果则说明有效粒子数 Meff不足,需要进行重采样,以解决粒子退化问题;所谓粒子退化问题是指,经过若干 次递推后,只有少数粒子的权重很大,而大部分粒子的权重非常接近于零的现象;由于权 重接近于零的粒子对状态的估计的贡献很小,所以当粒子退化严重时,有很多计算和存储 资源都被浪费了。其它步骤及参数与具体实施方式一至六之一相同。
采用以下实施例验证本发明的有益效果:
实施例一
步骤一、完成两个变量的初始化,首先将确定用户的初始状态,在用户初始状态的 领域内均匀分布10000个初始状态向量粒子均匀分布在x0领域 内;同时设置迭代次数的计数值为k=1;根据用户初始状态向量粒子可以得到k=1 时刻用户状态向量粒子
其中,
步骤二、当k≠1时,根据k-1时刻的用户状态向量粒子更新m个在第k时刻状 态向量粒子
是粒子滤波的动态处理模型,以适应GNSS接收机的运动;其中,状态向量粒子 也就是用户在k时隙的状态,也是最终的定位结果;所有坐 标都是在ECEF(Earth-Centered and Earth-Fixed,地心固定坐标系)坐标系下给出, (xU,k,yU,k,zU,k)代表用户在k个时刻状态向量粒子的三维位置;δtU,k用户在k个 时刻状态与第n颗卫星的时钟偏差;模型输入为:分别表示xU,k,yU,k,zU,k,δtU,k在单位时隙中的变化量,这里的单位时 隙为1秒;表示用户在第k时刻的状态向量粒子的过程噪声,服从均值为零的 高斯分布;
步骤三、根据步骤二计算的计算出用户与第n颗卫星之间的观测伪距由 下式给出:
是接收机和卫星之间的观测伪距,(xn,k,yn,k,zn,k)(n=1,2,...,N)代表第 n颗卫星在第k个时刻的位置坐标,c表示光速;
(1)通过用户在第k时刻的状态向量粒子计算得到观测伪距
其中,是接收机和卫星之间的观测伪距,(xn,k,yn,k,zn,k)(n=1,2,...,N)代 表第n颗卫星在第k个时刻的位置坐标,c表示光速;(xU,k,yU,k,zU,k,δtU,k)代表用 户在第k时刻的状态向量粒子的状态,包括用户在第k时刻的状态向量粒子三 维位置(xU,k,yU,k,zU,k)和相对于第n颗卫星的时钟偏差δtU,k;
(2)真实伪距ρn,k由接收机的伪距测量电路直接测量得到,假设可见卫星的数目为N, 信号捕获和跟踪的影响被忽略;在第k个时刻,测得用户和第n颗卫星之间的观测伪距与真实伪距ρn,k的关系为:
εn,k表示在k时刻接收机测量与第n颗卫星的伪距时产生的测量误差;
步骤四、根据步骤三得到的观测伪距更新观测向量
其中,用户与第n颗卫星的伪距观测误差为ρn,k表示在第k个时 刻,第n颗卫星与用户之间的真实伪距,N为可见卫星的个数;
步骤五、根据更新观测向量
计算未归一化的粒子权重
其中,M为在第k时刻的状态向量粒子数,取10000,N为可见卫星的个数,取5,σ 为的标准差,σ=5.9;
(1)建立离散时间的粒子滤波模型为:
观测向量为
(2)计算未归一化的粒子权重的过程为:
其中,未归一化的粒子权重可由下式得到:
上式中,q(·)为重要密度函数(这个是为了推导方便,人为定义的一个函数,没有明 确的意义),每一个状态向量粒子都来自对q(·)的完备采样,为已知 时的概率密度函数,为已知时的概率密度函数;为了 得到方差最小的未归一化粒子权重,应该选取的后验概率密度 即已知和时的概率密度,但在大多数实际情况中是 无法实现的,可以采用近似表示,得到了未归一化粒子权重的表达式:
其中,经过重采样后,我们可以认为
M=10000不失一般性地假设第N颗卫星是欺骗卫星,从而得到了未归一化粒子权重 的最终表达:
其中,M为在第k时刻的状态向量粒子数,取10000,N为可见卫星的个数,取5,σ 为的标准差,σ=5.9;
归一化后为:
步骤六、根据未归一化的粒子权重计算观测向量的模值的最小值θ:
假设有多于四颗可见的真实卫星,还有一颗欺骗卫星;这颗欺骗卫星会给接收机提供 一个自身坐标(xf,k,yf,k,zf,k)和假的伪距从而提出了对欺骗的侦测和抑制方法, 由下式给出:
ρf,k是接收机和卫星之间的真实伪距,是欺骗卫星提供的假伪距,△ρf,k是二者的偏 差;一般而言,定位误差与|△ρf,k|正相关;(本公式是理论上的推导,说明本方法的合理 性)
前面已经提到,就是测量伪距与理论计算得到的伪距之间的差值;在没有欺骗卫 星的条件下,该向量中的每一个数都不会很大;如果存在欺骗卫星,对应的伪距误差就会 很大,从而导致的二范数很大,也就导致很大;所以,通过监测来侦测欺 骗卫星,具体如下:
这里的M为粒子数,取10000,根据定义监测变量θ为:
当存在欺骗卫星时且可见卫星数为N=5时,θ服从非中心卡方分布,其概率密度函 数为:
在N>5的条件下,gf(θ)可以近似为:
IT(·)为T阶修正的贝塞尔函数,
如果不存在欺骗卫星,则θ服从中心卡方分布,自由度为可见卫星数减4,即N-4:
其中,Γ(·)是gamma函数,λ为卡方分布的非中心度参数,其取值与可见卫星个数N 有关,也与伪距有关;由下面的方法得出:
bk=[0 0 ... bs,k ... 0]T
在第k个时刻,这里1n,k是当前接收机指向第n颗卫星的单位方向矢量,I表示单位 矩阵,bs,k为第n颗卫星的附加伪距,bk为对应的附加伪距矩阵,Gk表示接收机与所有 可见卫星的单位方向矩阵,T代表矩阵的转置,是Gk的转置;
可见,存在欺骗卫星的情况下和不存在欺骗卫星的情况下,θ服从不同的分布,可以 通过该分布来判断是否存在欺骗卫星;因此,可以设置一个门限γ,当θ>γ时认为存在 欺骗卫星,反之则认为没有;这里的γ与可见卫星数N有关,N=5时,γ取75,其他情 况下γ取值可用试验方法得到;依据T阶修正的贝塞尔函数IT(·)说明本方法的合理性如 图3所示;
步骤七、判断θ与γ的大小,如果θ>γ,根据修正这 里的γ取75,得到:
其中,b为误差下界取6,d为误差上界取8,为观测误差,Pf(·)为误差抑制函 数;
然后根据修正后的再重新计算未归一化的粒子权重使得θ≤γ;
当不存在欺骗卫星时,该方法有一定的概率发生误判,这一误判的概率被称为虚警概 率Pf,与门限γ有关;虚警概率Pf与门限γ的关系由下面的式子给出:
相应的,发现概率为:
其中,g(θ):不存在欺骗卫星的情况下,监测变量θ的概率密度函数;PD:发现欺 骗卫星的概率;gf(θ):存在欺骗卫星的情况下,监测变量θ的概率密度函数;
在侦测到欺骗卫星后,可以通过基于稳健估计的改进方法来对其进行抑制;稳健估计 的目的在于,在一些普遍假设的条件下,给出最贴近真实值的估计;传统的估计方法,比 如最小二乘法,认为每一个观测数据都服从先验条件,所以,某一个错误的观测值会带来 很大的误差;而稳健估计则不同,它会去除某些异常观测值,从而消除其影响;最大似然 类型的估计(M-estimate),高阶统计的线性组合(L-estimate),等,是基本的稳健估计方法; 其中,最大似然比估计总是被用在参数估计上;
本发明基于传统的稳健估计方法,提出了改进,具体处理方法如下:
首先定义误差抑制函数Pf(·):
其中,b为误差下界,d为误差上界,e为观测误差,在计算得到未归一化粒子权重后,对观测误差进行修正,得到方法如下:
经过这一处理后:
(1)若误差值的绝对值超过d(可以取d=8,b=6),则说明对应的卫星是一颗假 卫星,将它修正为0,即它就不会对后续计算产生 任何影响了,也就是说不会影响到定位结果,所以假卫星被完全抑制,
(2)若的绝对值在b和d之间,则修正后减小为b,即
(3)若的绝对值小于b,则其大小则保持不变,即 然后再重新计算未归一化粒子权重,得到归一化粒 子权重这一部分的作用在于,抑制了欺骗卫星对定位结果的干扰;如果θ≤γ,则 直接跳到步骤八;
步骤八、根据未归一化的粒子权重计算k时刻定位结果
考虑最小均方估计,则对于的估计为
其中,为归一化粒子权重;即作为第k个时刻的定位结果输出;
步骤九、将步骤二中共得到M个粒子在k时刻的状态向量粒子由于存在粒子退 化现象,部分状态粒子将失效,剩余的有效状态粒子的数目Meff,根据步骤五中的未归 一化的粒子权重计算得到剩余的有效向量粒子的数目Meff:
如果则转到步骤十一,如果跳到步骤十,进行重采样;
步骤十、如果采取M个在k时刻的状态向量粒子对集合重采样,采到的概率为归一化粒子权重重新采到的归一化粒子权重直至使得 为止的具体过程为:
对集合中的元素进行均匀采样,依次抽取M次,得到k时刻M个新的状 态向量粒子,用k时刻M个新的状态向量粒子来取代之前的k时刻的M个状态向量粒子 其中M个状态向量粒子抽取的次数为M次,如果则说明有效 粒子数Meff不足,需要进行重采样,以解决粒子退化问题;所谓粒子退化问题是指,经 过若干次递推后,只有少数粒子的权重很大,而大部分粒子的权重非常接近于零的现象; 由于权重接近于零的粒子对状态的估计的贡献很小,所以当粒子退化严重时,有很多计算 和存储资源都被浪费了;
步骤十一、如果利用步骤二中的概率为的在k时刻状态向量粒子 将k时刻加1后得到用户在k+1时刻的用户状态向量粒子跳到步骤二,开 始下一个时刻的定位过程;
仿真结果表明,该发明可以有效侦测和抑制欺骗卫星的干扰,首先给出侦测效果的仿真 结果,仿真参数设置如下表:
实施例二
步骤一、完成两个变量的初始化,首先将确定用户的初始状态,在用户初始状态的 领域内均匀分布10000个初始状态向量粒子均匀分布在x0领域 内;同时设置迭代次数的计数值为k=1;根据用户初始状态向量粒子可以得到k=1 时刻用户状态向量粒子
其中,
步骤二、当k≠1时,根据k-1时刻的用户状态向量粒子更新m个在第k时刻状 态向量粒子
是粒子滤波的动态处理模型,以适应GNSS接收机的运动;其中,状态向量粒子 也就是用户在k时隙的状态,也是最终的定位结果;所有坐 标都是在ECEF(Earth-Centered and Earth-Fixed,地心固定坐标系)坐标系下给出, (xU,k,yU,k,zU,k)代表用户在k个时刻状态向量粒子的三维位置;δtU,k用户在k个 时刻状态与第n颗卫星的时钟偏差;模型输入为:分别表示xU,k,yU,k,zU,k,δtU,k在单位时隙中的变化量,这里的单位时 隙为1秒;表示用户在第k时刻的状态向量粒子的过程噪声,服从均值为零的 高斯分布;
步骤三、根据步骤二计算的计算出用户与第n颗卫星之间的观测伪距由 下式给出:
是接收机和卫星之间的观测伪距,(xn,k,yn,k,zn,k)(n=1,2,...,N)代表第 n颗卫星在第k个时刻的位置坐标,c表示光速;
(1)通过用户在第k时刻的状态向量粒子计算得到观测伪距
其中,是接收机和卫星之间的观测伪距,(xn,k,yn,k,zn,k)(n=1,2,...,N)代 表第n颗卫星在第k个时刻的位置坐标,c表示光速;(xU,k,yU,k,zU,k,δtU,k)代表用 户在第k时刻的状态向量粒子的状态,包括用户在第k时刻的状态向量粒子三 维位置(xU,k,yU,k,zU,k)和相对于第n颗卫星的时钟偏差δtU,k;
(2)真实伪距ρn,k由接收机的伪距测量电路直接测量得到,假设可见卫星的数目为N, 信号捕获和跟踪的影响被忽略;在第k个时刻,测得用户和第n颗卫星之间的观测伪距与真实伪距ρn,k的关系为:
εn,k表示在k时刻接收机测量与第n颗卫星的伪距时产生的测量误差;
步骤四、根据步骤三得到的观测伪距更新观测向量
其中,用户与第n颗卫星的伪距观测误差为ρn,k表示在第k个时 刻,第n颗卫星与用户之间的真实伪距,N为可见卫星的个数;
步骤五、根据更新观测向量
计算未归一化的粒子权重
其中,M为在第k时刻的状态向量粒子数,取10000,N为可见卫星的个数,取5,σ 为的标准差,σ=5.9;
(1)建立离散时间的粒子滤波模型为:
观测向量为代表了测量得到的伪距与理论计 算得到的伪距之间的差值;ε(m)k表示第k时刻的状态向量粒子的观测噪声;h(m)(·) 代表与yk(m)之间的函数关系;
(2)计算未归一化的粒子权重的过程为:
其中,未归一化的粒子权重可由下式得到:
上式中,q(·)为重要密度函数(这个是为了推导方便,人为定义的一个函数,没有明 确的意义),每一个状态向量粒子都来自对q(·)的完备采样,为已知 时的概率密度函数,为已知时的概率密度函数;为了 得到方差最小的未归一化粒子权重,应该选取的后验概率密度 即已知和时的概率密度,但在大多数实际情况中是 无法实现的,可以采用近似表示,得到了未归一化粒子权重的表达式:
其中,经过重采样后,我们可以认为
M=10000不失一般性地假设第N颗卫星是欺骗卫星,从而得到了未归一化粒子权重 的最终表达:
其中,M为在第k时刻的状态向量粒子数,取10000,N为可见卫星的个数,取8,σ 为的标准差,σ=5.9;
归一化后为:
步骤六、根据未归一化的粒子权重计算观测向量的模值的最小值θ:
假设有多于四颗可见的真实卫星,还有一颗欺骗卫星;这颗欺骗卫星会给接收机提供 一个自身坐标(xf,k,yf,k,zf,k)和假的伪距从而提出了对欺骗的侦测和抑制方法, 由下式给出:
ρf,k是接收机和卫星之间的真实伪距,是欺骗卫星提供的假伪距,△ρf,k是二者的偏 差;一般而言,定位误差与|△ρf,k|正相关;(本公式是理论上的推导,说明本方法的合理 性)
前面已经提到,就是测量伪距与理论计算得到的伪距之间的差值;在没有欺骗卫 星的条件下,该向量中的每一个数都不会很大;如果存在欺骗卫星,对应的伪距误差就会 很大,从而导致的二范数很大,也就导致很大;所以,通过监测来侦测欺 骗卫星,具体如下:
这里的M为粒子数,取10000,根据定义监测变量θ为:
当存在欺骗卫星时且可见卫星数为N=5时,θ服从非中心卡方分布,其概率密度函数 为:
在N>5的条件下,gf(θ)可以近似为:
IT(·)为T阶修正的贝塞尔函数,
如果不存在欺骗卫星,则θ服从中心卡方分布,自由度为可见卫星数减4,即N-4:
其中,Γ(·)是gamma函数,λ为卡方分布的非中心度参数,其取值与可见卫星个数N 有关,也与伪距有关;由下面的方法得出:
bk=[0 0 ... bs,k ... 0]T
在第k个时刻,这里1n,k是当前接收机指向第n颗卫星的单位方向矢量,I表示单位 矩阵,bs,k为第n颗卫星的附加伪距,bk为对应的附加伪距矩阵,Gk表示接收机与所有 可见卫星的单位方向矩阵,T代表矩阵的转置,是Gk的转置;
可见,存在欺骗卫星的情况下和不存在欺骗卫星的情况下,θ服从不同的分布,可以 通过该分布来判断是否存在欺骗卫星;因此,可以设置一个门限γ,当θ>γ时认为存在 欺骗卫星,反之则认为没有;这里的γ与可见卫星数N有关,N=8时,γ取250,其他情 况下γ取值可用试验方法得到;依据T阶修正的贝塞尔函数IT(·)说明本方法的合理性如 图4所示;
步骤七、判断θ与γ的大小,如果θ>γ,根据修正这 里的γ取250,得到:
其中,b为误差下界取6,d为误差上界取8,为观测误差,Pf(·)为误差抑制函 数;
然后根据修正后的再重新计算未归一化的粒子权重使得θ≤γ;
当不存在欺骗卫星时,该方法有一定的概率发生误判,这一误判的概率被称为虚警概 率Pf,与门限γ有关;虚警概率Pf与门限γ的关系由下面的式子给出:
相应的,发现概率为:
其中,g(θ):不存在欺骗卫星的情况下,监测变量θ的概率密度函数;PD:发现欺 骗卫星的概率;gf(θ):存在欺骗卫星的情况下,监测变量θ的概率密度函数;
在侦测到欺骗卫星后,可以通过基于稳健估计的改进方法来对其进行抑制;稳健估计 的目的在于,在一些普遍假设的条件下,给出最贴近真实值的估计;传统的估计方法,比 如最小二乘法,认为每一个观测数据都服从先验条件,所以,某一个错误的观测值会带来 很大的误差;而稳健估计则不同,它会去除某些异常观测值,从而消除其影响;最大似然 类型的估计(M-estimate),高阶统计的线性组合(L-estimate),等,是基本的稳健估计方法; 其中,最大似然比估计总是被用在参数估计上;
本发明基于传统的稳健估计方法,提出了改进,具体处理方法如下:
首先定义误差抑制函数Pf(·):
其中,b为误差下界,d为误差上界,e为观测误差,在计算得到未归一化粒子权重后,对观测误差进行修正,得到方法如下:
经过这一处理后:
(1)若误差值的绝对值超过d(可以取d=8,b=6),则说明对应的卫星是一颗假 卫星,将它修正为0,即它就不会对后续计算产生 任何影响了,也就是说不会影响到定位结果,所以假卫星被完全抑制,
(2)若的绝对值在b和d之间,则修正后减小为b,即
(3)若的绝对值小于b,则其大小则保持不变,即 然后再重新计算未归一化粒子权重,得到归一化粒 子权重这一部分的作用在于,抑制了欺骗卫星对定位结果的干扰;如果θ≤γ,则 直接跳到步骤八;
步骤八、根据未归一化的粒子权重计算k时刻定位结果
考虑最小均方估计,则对于的估计为
其中,为归一化粒子权重;即作为第k个时刻的定位结果输出;
步骤九、将步骤二中共得到M个粒子在k时刻的状态向量粒子由于存在粒子退 化现象,部分状态粒子将失效,剩余的有效状态粒子的数目Meff,根据步骤五中的未归 一化的粒子权重计算得到剩余的有效向量粒子的数目Meff:
如果则转到步骤十一,如果跳到步骤十,进行重采样;
步骤十、如果采取M个在k时刻的状态向量粒子对集合重采样,采到的概率为归一化粒子权重重新采到的归一化粒子权重直至使得 为止的具体过程为:
对集合中的元素进行均匀采样,依次抽取M次,得到k时刻M个新的状 态向量粒子,用k时刻M个新的状态向量粒子来取代之前的k时刻的M个状态向量粒子 其中M个状态向量粒子抽取的次数为M次;如果则说明有效 粒子数Meff不足,需要进行重采样,以解决粒子退化问题;所谓粒子退化问题是指,经 过若干次递推后,只有少数粒子的权重很大,而大部分粒子的权重非常接近于零的现象; 由于权重接近于零的粒子对状态的估计的贡献很小,所以当粒子退化严重时,有很多计算 和存储资源都被浪费了;
步骤十一、如果利用步骤二中的概率为的在k时刻状态向量粒子 将k时刻加1后得到用户在k+1时刻的用户状态向量粒子跳到步骤二,开 始下一个时刻的定位过程;
仿真结果表明,该发明可以有效侦测和抑制欺骗卫星的干扰,首先给出侦测效果的仿 真结果,仿真参数设置如下表:
机译: IGBT一种基于辅助粒子滤波的IGBT剩余使用寿命的估算方法
机译: 一种基于粒子滤波的鲁棒目标跟踪方法及系统
机译: 一种基于粒子滤波框架的鲁棒目标跟踪方法及系统