CN103983986A - 一种基于粒子滤波的改进型raim抗欺骗式干扰方法 - Google Patents

Abstract

一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，本发明涉及改进型RAIM抗欺骗式干扰方法。本发明是要解决针对单颗卫星的测量失效，控制解算流程误导接收机定位以及RAIM忽略了残余矢量之间的相关性和相似性的问题，而提出的一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法。该方法是通过1、组成2、得到ρ_n，和卫星的三维坐标；3、计算4、计算y^(m)；5、计算w^(m)；6、选取w^max；7计算8、计算最大可见星SL^max；9、计算判决门限γ；10、判断存在欺骗卫星；11、估计干扰卫星的序号；12、若欺骗卫星标志F＝1，去掉干扰，进行定位解算，若F≠1进行定位解算等步骤实现的。本发明应用于改进型RAIM抗欺骗式干扰领域。

Description

一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法

技术领域

本发明涉及基于粒子滤波的改进型RAIM抗欺骗式干扰方法。

背景技术

定位服务是GNSS(全球导航卫星系统)所提供的最基本服务，同时它也是欺骗式干扰攻击的首选目标。干扰机通过转发或者伪造出带有欺骗星历信息和伪距信息的虚假信号，然后将这些信号传输到接收机空间中相应的频段上。通常，干扰方会利用压制式干扰破坏真实GNSS信号链路，然后增大干扰信号的功率以提高其被捕获的概率。被接收机捕获的欺骗式干扰信号可以破坏甚至控制解算流程，进而成功的误导接收机定位。

通过观察发现，有意的欺骗式干扰攻击和无意的GNSS接收机伪距测量失效非常相似。伪距测量失效是指GNSS接收机在测量可见星伪距时，出现稳定的测量偏差，进而破坏定位解算过程。针对伪距测量失效的情况，常用RAIM(收机自体完好性监控)实现失效卫星的检测和识别。

传统RAIM技术最早是基于最小二乘解算算法，通过二乘法的解算残余实现GNSS完整性检测。因为GNSS系统中卫星失效的概率很低，双星甚至多颗卫星同时失效的概率几乎为零，因而传统RAIM都是针对单星失效的场景进行讨论。它可以在可见星数量大于等于5颗时，实现卫星失效情况的检测，在可见卫星数量大于等于6颗时，实现失效卫星的序号的识别。类似的，传统RAIM技术也可以在可见星数量大于等于5和6的时候，实现干扰的检测和干扰卫星的识别。

在二十年前，很少有人考虑多颗干扰卫星的情况，因为当时技术限制了这一情况的出现。但随着集成电子技术和传感器技术的迅速发展，一个综合的干扰系统可以准确测量GNSS接收机的位置，并且感知其所有的可见星信息。那么干扰机可以根据当前的已知信息，模拟产生或者转发一个或者多个虚假的GNSS信号，因而提出抗多星欺骗式干扰的方案非常必要。传统RAIM最初只是针对单颗卫星的测量失效，尽管有人针对多颗卫星测量失效的情况对RAIM技术进行修正，但是其忽略了残余矢量之间的相关性和相似性，这就给干扰源以可乘之机。干扰机可以利用这一漏洞成功的避开干扰检测或者扰乱干扰识别过程。本发明正是为了弥补RAIM方法的这些缺陷。

发明内容

本发明的目的是为了解决接收机捕获的欺骗式干扰信号破坏导致控制解算流程误导接收机定位、传统RAIM最初只是针对单颗卫星的测量失效，以及RAIM忽略了残余矢量之间的相关性和相似性，而提出了一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法。

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

步骤一、接收机状态的估计值为表示接收机三维坐标的估计值，表示接收机时钟偏差的估计值，在的邻域内得到M个粒子，组成集合其中x^(m)＝(x_U,m,y_U,m,z_U,m,δt_U,m)表示第m个接收机状态粒子，(x_U,m,y_U,m,z_U,m)表示接收机的三维坐标，δt_U,m表示接收机时钟偏差，m＝1，2....M；

步骤二、设当前可见卫星数为N，测量每一颗卫星与接收机的伪距，得到第n颗卫星与接收机的测量伪距ρ_n(n＝1,2...N)，以及卫星的三维坐标(x_n,y_n,z_n)；

步骤三、根据步骤一得到的第m个接收机状态粒子x^(m)，及步骤二得到的第n颗卫星的三维坐标(x_n,y_n,z_n)，得到第n颗卫星与第m个粒子的计算伪距

其中，c代表光速；

步骤四、根据步骤二得到的第n颗卫星与接收机的测量伪距ρ_n(n＝1,2...N)，及步骤三得到的第n颗卫星与第m个粒子的计算伪距计算第m个粒子的观测向量y^(m)＝[e_1,m,...,e_n,m,...,e_N,m]^T，其中e_n,m表示第m个粒子与第n颗卫星的伪距残差：

步骤五、根据步骤四得到的第m个粒子的观测向量y^(m)，计算第m个粒子的未归一化加权值w^(m)：

$w^{\left(m\right)}=\frac{1}{M}{\left(\frac{1}{\sqrt{2{\mathrm{\π\σ}}^2}}\right)}^N\exp [-\frac{\left|\right|y^{\left(m\right)}\left|\right|}{{2\mathrm{\σ}}^2}];$

其中，σ表示无噪声条件下e_n,m的标准差，通常取5.9；

步骤六、根据步骤五得到的第m个粒子的未归一化加权值w^(m)，选取最大的粒子加权值w^max：

w^max＝max(w^(m))，

最大的粒子加权值w^max分别与最大观测向量y^max，最大接收机状态粒子x^max，对应的最大粒子偏差d^max相对应，其中，

$d^{\max }=\hat{x}-x^{\max },$

其中表示接收机状态的估计值；

步骤七、根据步骤六得到的最大观测向量y^max，计算附加伪距的估计值

步骤八、根据步骤七得到的附加伪距的估计值计算最大可见星联合斜率SL^max：

步骤九、根据步骤八得到的最大可见星联合斜率SL^max，计算判决门限γ；

步骤十、根据步骤九得到的判决门限γ和步骤六得到的最大观测向量y^max，若‖y^max‖≥γ则说明存在欺骗卫星，欺骗卫星标志F＝1，若不存在则跳到步骤十二，欺骗卫星标志F＝0；

步骤十一、则根据步骤六得到的最大观测向量为y^max、步骤七得到的附加伪距的估计值和残余矢量矩阵S，估计干扰卫星的序号

步骤十二、若欺骗卫星标志F＝1，去掉步骤十一的估计干扰卫星的序号的干扰，然后进行定位解算，若欺骗卫星标志F≠1直接进行定位解算；即完成了一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法。

发明效果

本发明提出了一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，能够在不增加任何硬件资源、不改变信号体制的条件下实现欺骗卫星的检测和识别。从而确保了定位精度的准确，保护了国家和人民的利益。所有坐标都是在ECEF(Earth-CenteredandEarth-Fixed，地心固定坐标系)坐标系下给出。

图2比较了传统RAIM和改进型RAIM对于双星干扰(ModeA)的干扰检测概率和干扰识别概率。可以看出，本文提出的改进型RAIM可以提供更好的抗干扰效果。ModeA通过两个卫星附加伪距的相互配合来降低‖y^max‖，使传统RAIM无法有效地检测到干扰攻击的存在，因而其检测概率比改进型RAIM的检测概率低。同时，当定位误差大于450m，检测概率都为1时，改进型RAIM也可以提供更好的检测性能，这也证明了改进的有效性。

图3给出了双星干扰(ModeB)的检测概率P_I和错误将第五号卫星识别为干扰的概率P_MI。从图中可以看出，传统RAIM的误识别概率P_MI随着定位误差的增大而增加，这是因为接收机首先假设只有一颗干扰卫星，而(0.65s¹+s⁶)和之间的高相似性会使接收机将第五号卫星误识别干扰星。而本文提出的改进型RAIM会对识别出的疑似干扰卫星进行验证，从而有效地降低了误识别概率，增加正确识别概率。

图4给出了双星干扰(ModeC)条件下的干扰识别概率P_I和将第六号和第七号卫星误识别为干扰的概率P_MI正确识别概率P_I和误识别概率P_MI都应该为0.5左右，改进型RAIM的正确识别概率随着定位误差的增加而增加。

在实施例ModeA中都可以解释。当干扰卫星的数量大于1颗时，可以看到，图2中传统RAIM的曲线，效果很差，说明它对多颗干扰卫星的情况无能为力。同时，在ModeA中，通过两个卫星附加伪距的相互配合来降低‖y^max‖，即两颗干扰卫星的残余矢量相关性很强，在这种情况下，本发明仍旧在很大概率下检测到了干扰。

附图说明

图1是具体实施方式一提出的一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法流程图；

图2是实施例提出的双星干扰(ModeA)的干扰检测概率和识别概率示意图；

图3是实施例提出的双星干扰(ModeB)的干扰检测概率和识别概率示意图；

图4是实施例提出的双星干扰(ModeC)的干扰检测概率和识别概率示意图。

具体实施方式

具体实施方式一：本实施方式的一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，具体是按照以下步骤制备的：

步骤一、接收机状态的估计值为表示接收机三维坐标的估计值，表示接收机时钟偏差的估计值，在的邻域内得到M个粒子，组成了集合其中x^(m)＝(x_U,m,y_U,m,z_U,m,δt_U,m)表示第m个接收机状态粒子，(x_U,m,y_U,m,z_U,m)表示接收机的三维坐标，δt_U,m表示接收机时钟偏差，m＝1，2....M；

步骤二、设当前可见卫星数为N，测量每一颗卫星与接收机的伪距，得到第n颗卫星与接收机的测量伪距ρ_n(n＝1,2...N)，以及卫星的三维坐标(x_n,y_n,z_n)；

步骤三、根据步骤一得到的第m个接收机状态粒子x^(m)，及步骤二得到的第n颗卫星的三维坐标(x_n,y_n,z_n)，得到第n颗卫星与第m个粒子的计算伪距

其中，c代表光速；

步骤四、根据步骤二得到的第n颗卫星与接收机的测量伪距ρ_n(n＝1,2...N)，及步骤三得到的第n颗卫星与第m个粒子的计算伪距计算第m个粒子的观测向量y^(m)＝[e_1,m,...,e_n,m,...,e_N,m]^T，其中e_n,m表示第m个粒子与第n颗卫星的伪距残差：

步骤五、根据步骤四得到的第m个粒子的观测向量y^(m)，计算第m个粒子的未归一化加权值w^(m)：

$w^{\left(m\right)}=\frac{1}{M}{\left(\frac{1}{\sqrt{2{\mathrm{\π\σ}}^2}}\right)}^N\exp [-\frac{\left|\right|y^{\left(m\right)}\left|\right|}{{2\mathrm{\σ}}^2}];$

其中，σ表示无噪声条件下e_n,m的标准差，通常取5.9；

步骤六、根据步骤五得到的第m个粒子的未归一化加权值w^(m)，选取最大的粒子加权值w^max：

w^max＝max(w^(m))，

最大的粒子加权值w^max分别与最大观测向量y^max，最大接收机状态粒子x^max，对应的最大粒子偏差d^max相对应，其中，

$d^{\max }=\hat{x}-x^{\max },$

其中表示接收机状态的估计值；

步骤七、根据步骤六得到的最大观测向量y^max，计算附加伪距的估计值

步骤八、根据步骤七得到的附加伪距的估计值计算最大可见星联合斜率SL^max：

步骤九、根据步骤八得到的最大可见星联合斜率SL^max，计算判决门限γ；

步骤十、根据步骤九得到的判决门限γ和步骤六得到的最大观测向量y^max，判断是否存在欺骗卫星：若‖y^max‖≥γ则说明存在欺骗卫星，欺骗卫星标志F＝1，若不存在则跳到步骤十二，欺骗卫星标志F＝0；

步骤十一、则根据步骤六得到的最大观测向量为y^max、步骤七得到的附加伪距的估计值和残余矢量矩阵S，估计干扰卫星的序号

步骤十二、若欺骗卫星标志F＝1，去掉步骤十一的估计干扰卫星的序号的干扰，然后进行定位解算，若欺骗卫星标志F≠1直接进行定位解算如图1；即完成了一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法。

本实施方式效果

本实施方式提出了一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，能够在不增加任何硬件资源、不改变信号体制的条件下实现欺骗卫星的检测和识别。从而确保了定位精度的准确，保护了国家和人民的利益。所有坐标都是在ECEF(Earth-Centered andEarth-Fixed，地心固定坐标系)坐标系下给出。

图2比较了传统RAIM和改进型RAIM对于双星干扰(ModeA)的干扰检测概率和干扰识别概率。可以看出，本文提出的改进型RAIM可以提供更好的抗干扰效果。ModeA通过两个卫星附加伪距的相互配合来降低‖y^max‖，使传统RAIM无法有效地检测到干扰攻击的存在，因而其检测概率比改进型RAIM的检测概率低。同时，当定位误差大于450m，检测概率都为1时，改进型RAIM也可以提供更好的检测性能，这也证明了改进的有效性。

图3给出了双星干扰(ModeB)的检测概率P_I和错误将第五号卫星识别为干扰的概率P_MI。从图中可以看出，传统RAIM的误识别概率P_MI随着定位误差的增大而增加，这是因为接收机首先假设只有一颗干扰卫星，而(0.65s¹+s⁶)和之间的高相似性会使接收机将第五号卫星误识别干扰星。而本文提出的改进型RAIM会对识别出的疑似干扰卫星进行验证，从而有效地降低了误识别概率，增加正确识别概率。

图4给出了双星干扰(ModeC)条件下的干扰识别概率P_I和将第六号和第七号卫星误识别为干扰的概率P_MI正确识别概率P_I和误识别概率P_MI都应该为0.5左右，改进型RAIM的正确识别概率随着定位误差的增加而增加。

在实施例ModeA中都可以解释。当干扰卫星的数量大于1颗时，可以看到，图2中传统RAIM的曲线，效果很差，说明它对多颗干扰卫星的情况无能为力。同时，在ModeA中，通过两个卫星附加伪距的相互配合来降低‖y^max‖，即两颗干扰卫星的残余矢量相关性很强，在这种情况下，本实施方式仍旧在很大概率下检测到了干扰。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤七中根据步骤六得到的最大观测向量y^max，计算附加伪距的估计值为：

$\hat{b}={\left(S^{T}S\right)}^{-1}{S}^T{y}^{\max },$

其中S＝(I-G(G^TG)^-1G^T)，为残余矢量矩阵，G代表卫星单位方向矢量矩阵，由下式得到：

$G=\left[\begin{array}{cc}{1}_1& 1\\ 1_2& 1\\ .& .\\ .& .\\ .& .\\ 1_N& 1\end{array}\right],$

1_n(n＝1,2,...,N)是地心指向第n颗卫星的单位方向矢量，I为4×4的单位矩阵。其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤八中根据步骤七得到的附加伪距的估计值计算最大可见星联合斜率SL^max为：

${\mathrm{SL}}^{\max }=\underset{n_1,n_2,...,n_{\mathrm{Ns}}}{\max }\left[\frac{\left|\right|{\left(d^{\max }\right)}_{\mathrm{sub}}^{\mathrm{1,2,3}}\left|\right|}{\left|\right|S_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}{\hat{b}}_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}\left|\right|}\right],$

其中，1≤n₁<n₂<...<n_Ns≤N，N代表可见卫星总数，Ns表示干扰卫星总数，n₁,n₂,...,n_Ns均为整数代表可见卫星的编号，表示d^max的前三个元素构成的子矩阵，表示S的n₁,n₂,...,n_Ns列构成的子矩阵，表示的n₁,n₂,...,n_Ns行构成的子矩阵。其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤九中根据步骤八得到的最大可见星联合斜率SL^max，计算判决门限γ为：

$\mathrm{\&gamma;}=\frac{\mathrm{PL}}{{\mathrm{SL}}^{\max }},$

其中PL为定位误差保护门限。其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤十一中根据步骤六得到的最大观测向量为y^max、步骤七得到的附加伪距的估计值和残余矢量矩阵S，估计干扰卫星的序号为：

$({\hat{s}}_1,{\hat{s}}_2,...{\hat{s}}_{\mathrm{Ns}})=\arg \underset{n_1,n_2,...,n_{\mathrm{Ns}}}{\max }\left[\frac{\left|\right|y^{\max }\&CenterDot;\left(S_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}{\hat{b}}_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}\right)\left|\right|}{\left|\right|y^{\max }\left|\right|\left|\right|\left(S_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}{\hat{b}}_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}\right)\left|\right|}\right].$ 其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤十二中定位解算具体过程为：

去掉欺骗卫星后，假设剩余的可用卫星数目为M(卫星数目为M为可见卫星N减去干扰卫星的个数)，设它们的编号为：s₁,s₂......,s_M(1≤1s₁<s₂<...<s_M≤N)，接收机的坐标为(x_u,y_u,z_u)，编号为s₁的卫星与接收机之间的测量伪距为ρs₁，则对该卫星存在如下的伪距方程：

${\mathrm{\&rho;}}_{s_1}=\sqrt{{(x_u-x_{s_1})}^2+{(y_u-y_{s_1})}^2+{(z_u-z_{s_1})}^2}+c\&times;(\mathrm{\&Delta;T}-{\mathrm{\&Delta;t}}_{s_1})$

其中，代表编号为s₁的卫星的三维坐标，代表该卫星时钟相对于GPS时系的偏差，c代表光速；这三者对于接收机来说都是已知量；ΔT为接收机时钟相对GPS时系的偏差，ΔT取正值表示接收机时钟相对于GPS时系超前，ΔT不取正值表示滞后；令b_u＝c×ΔT，则

${\mathrm{\&rho;}}_{s_1}=\sqrt{{(x_u-x_{s_1})}^2+{(y_u-y_{s_1})}^2+{(z_u-z_{s_1})}^2}+b_u-c{\&times;\mathrm{\&Delta;t}}_{s_1}$

式中共有四个未知数，即接收机的三维坐标(x_u,y_u,z_u)和时钟偏差造成的伪距偏差b_u；对另外M-1颗卫星列伪距方程，即得到M个伪距方程，求解出接收机的三维坐标(x_u,y_u,z_u)和时钟偏差造成的伪距偏差b_u，得到接收机状态；其中，第M个方程是 ${\mathrm{\&rho;}}_{s_M}=\sqrt{{(x_u-x_{\mathrm{sM}})}^2+{(y_u-y_{\mathrm{sM}})}^2+{(z_u-z_{\mathrm{sM}})}^2}+b_u-c{\&times;\mathrm{\&Delta;t}}_{s_M}.$ 其它步骤及参数与具体实施方式一至五之一相同。

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

实施例一：

本实施例一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，具体是按照以下步骤制备的：

步骤一、接收机状态的估计值为表示接收机三维坐标的估计值，在后续仿真中设为(-2654200，3564100，4565300，0)表示接收机时钟偏差的估计值，在的邻域内得到M个粒子，组成了集合其中x^(m)＝(x_U,m,y_U,m,z_U,m,δt_U,m)表示第m个接收机状态粒子，(x_U,m,y_U,m,z_U,m)表示接收机的三维坐标，δt_U,m表示接收机时钟偏差，m＝1，2....M，M取10000；

步骤二、设当前可见卫星数为N，测量每一颗卫星与接收机的伪距，得到第n颗卫星与接收机的测量伪距ρn(n＝1,2...N)，以及该卫星的三维坐标(x_n,y_n,z_n)；假设当前接收机的可见星为8颗，可见卫星的位置的三维坐标(x_n,y_n,z_n)如表1所示：

表1可见星的位置坐标

卫星序号	x(m)	y(m)	z(m)
				1	-16763868	-4481176	20106830
2	-19725161	16121107	7280333
				3	-19799420	5599382	16634207
4	6523912	25435850	2950164
				5	-26479541	-3394195	574679
6	5506	15323904	21907827
				7	-15310052	21353100	2595262
8	-13237017	5964854	11241253

步骤三、根据步骤一得到的第m个接收机状态粒子x^(m)，及步骤二得到的第n颗卫星的三维坐标(x_n,y_n,z_n)，得到第n颗卫星与第m个粒子的计算伪距

其中，c代表光速；

步骤四、根据步骤二得到的第n颗卫星与接收机的测量伪距ρ_n(n＝1,2...N)，及步骤三得到的第n颗卫星与第m个粒子的计算伪距计算第m个粒子的观测向量y^(m)＝[e_1,m,...,e_n,m,...,e_N,m]^T，其中e_n,m表示第m个粒子与第n颗卫星的伪距残差：

步骤五、根据步骤四得到的第m个粒子的观测向量y^(m)，计算第m个粒子的未归一化加权值w^(m)：

$w^{\left(m\right)}=\frac{1}{M}{\left(\frac{1}{\sqrt{2{\mathrm{\&pi;\&sigma;}}^2}}\right)}^N\exp [-\frac{\left|\right|y^{\left(m\right)}\left|\right|}{{2\mathrm{\&sigma;}}^2}];$

其中，σ表示无噪声条件下e_n,m的标准差，通常取5.9；

表2仿真参数

步骤六、根据步骤五得到的第m个粒子的未归一化加权值w^(m)，选取最大的粒子加权值w^max：

w^max＝max(w^(m))，

最大的粒子加权值w^max分别与最大观测向量y^max，最大接收机状态粒子x^max，对应的最大粒子偏差d^max相对应，其中，

$d^{\max }=\hat{x}-x^{\max },$

其中表示接收机状态的估计值；

步骤七、根据步骤六得到的最大观测向量y^max，计算附加伪距的估计值

$\hat{b}={\left(S^{T}S\right)}^{-1}{S}^T{y}^{\max },$

其中S＝(I-G(G^TG)^-1G^T)，为残余矢量矩阵，G代表卫星单位方向矢量矩阵，由下式得到：

$G=\left[\begin{array}{cc}{1}_1& 1\\ 1_2& 1\\ .& .\\ .& .\\ .& .\\ 1_N& 1\end{array}\right],$

1_n(n＝1,2,...,N)是地心指向第n颗卫星的单位方向矢量，I为4×4的单位矩阵；

步骤八、根据步骤七得到的附加伪距的估计值计算最大可见星联合斜率SL^max：

${\mathrm{SL}}^{\max }=\underset{n_1,n_2,...,n_{\mathrm{Ns}}}{\max }\left[\frac{\left|\right|{\left(d^{\max }\right)}_{\mathrm{sub}}^{\mathrm{1,2,3}}\left|\right|}{\left|\right|S_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}{\hat{b}}_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}\left|\right|}\right],$

其中，1≤n₁<n₂<...<n_Ns≤N，N代表可见卫星总数，Ns表示干扰卫星总数，n₁,n₂,...,n_Ns均为整数代表可见卫星的编号，表示d^max的前三个元素构成的子矩阵，表示S的n₁,n₂,...,n_Ns列构成的子矩阵，表示的n₁,n₂,...,n_Ns行构成的子矩阵

步骤九、根据步骤八得到的最大可见星联合斜率SL^max，计算判决门限γ：

$\mathrm{\&gamma;}=\frac{\mathrm{PL}}{{\mathrm{SL}}^{\max }},$

其中PL为定位误差保护门限，设定γ为100米；

步骤十、根据步骤九得到的判决门限γ和步骤六得到的最大观测向量y^max，判断是否存在欺骗卫星：若‖y^max‖≥γ则说明存在欺骗卫星，欺骗卫星标志F＝1，若不存在则跳到步骤十二，欺骗卫星标志F＝0；

步骤十一、则根据步骤六得到的最大观测向量为y^max、步骤七得到的附加伪距的估计值和残余矢量矩阵S，估计干扰卫星的序号为：

$({\hat{s}}_1,{\hat{s}}_2,...{\hat{s}}_{\mathrm{Ns}})=\arg \underset{n_1,n_2,...,n_{\mathrm{Ns}}}{\max }\left[\frac{\left|\right|y^{\max }\&CenterDot;\left(S_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}{\hat{b}}_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}\right)\left|\right|}{\left|\right|y^{\max }\left|\right|\left|\right|\left(S_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}{\hat{b}}_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}\right)\left|\right|}\right];$

步骤十二、若欺骗卫星标志F＝1，去掉步骤十一的估计干扰卫星的序号的干扰，然后进行定位解算，若欺骗卫星标志F≠1直接进行定位解算；

去掉欺骗卫星后，假设剩余的可用卫星数目为M，设它们的编号为：s₁,s₂......,s_M(1≤1s₁<s₂<...<s_M≤N)，接收机的坐标为(x_u,y_u,z_u)，编号为s₁的卫星与接收机之间的测量伪距为则对该卫星存在如下的伪距方程：

${\mathrm{\&rho;}}_{s_1}=\sqrt{{(x_u-x_{s_1})}^2+{(y_u-y_{s_1})}^2+{(z_u-z_{s_1})}^2}+c\&times;(\mathrm{\&Delta;T}-{\mathrm{\&Delta;t}}_{s_1})$

其中，代表编号为s₁的卫星的三维坐标，代表该卫星时钟相对于GPS时系的偏差，c代表光速；这三者对于接收机来说都是已知量；ΔT为接收机时钟相对GPS时系的偏差，ΔT取正值表示接收机时钟相对于GPS时系超前，ΔT不取正值表示滞后；令b_u＝c×ΔT，则

${\mathrm{\&rho;}}_{s_1}=\sqrt{{(x_u-x_{s_1})}^2+{(y_u-y_{s_1})}^2+{(z_u-z_{s_1})}^2}+b_u-c{\&times;\mathrm{\&Delta;t}}_{s_1}$

式中共有四个未知数，即接收机的三维坐标(x_u,y_u,z_u)和时钟偏差造成的伪距偏差b_u；对另外M-1颗卫星列伪距方程，即得到M个伪距方程，求解出接收机的三维坐标(x_u,y_u,z_u)和时钟偏差造成的伪距偏差b_u即接收机状态；其中，第M个方程是 ${\mathrm{\&rho;}}_{s_M}=\sqrt{{(x_u-x_{\mathrm{sM}})}^2+{(y_u-y_{\mathrm{sM}})}^2+{(z_u-z_{\mathrm{sM}})}^2}+b_u-c{\&times;\mathrm{\&Delta;t}}_{s_M};$

假设干扰卫星有两颗，分三种干扰模式进行仿真，干扰卫星的编号和附加伪距见表3：

表3双星干扰攻击的仿真参数

这里的干扰模式可用下面的公式表示：

$\left\{\begin{array}{cc}\min \left|\right|s^{s_1}{b}^{s_1}+s^{s_2}{b}^{s_2}\left|\right|& \mathrm{ModeA}\\ \max \left|\mathrm{COV}(s^{s_1}{b}^{s_1}+s^{s_2}{b}^{s_2},s^{\mathrm{ms}}{b}^{\mathrm{ms}})\right|& \mathrm{ModeB}\\ \max \left|\mathrm{COV}(s^{s_1}{b}^{s_1}+s^{s_2}{b}^{s_2},s^{{\mathrm{ms}}_1}{b}^{{\mathrm{ms}}_1}+s^{{\mathrm{ms}}_2}{b}^{{\mathrm{ms}}_2})\right|& \mathrm{ModeC}\end{array}\right.$

ModeA旨在避开干扰检测，而ModeB和ModeC则希望扰乱干扰识别的过程；在ModeB中，干扰机精确控制和使y^max与s^msb^ms相似，进而误导接收将第ms号卫星识别为干扰；类似地，在ModeC中，干扰机误导接收将第ms₁和ms₂号卫星识别为干扰；

图2比较了传统RAIM和改进型RAIM对于双星干扰(ModeA)的干扰检测概率和干扰识别概率；可以看出，本文提出的改进型RAIM可以提供更好的抗干扰效果；ModeA通过两个卫星附加伪距的相互配合来降低‖y^max‖，使传统RAIM无法有效地检测到干扰攻击的存在，因而其检测概率比改进型RAIM的检测概率低；同时，当定位误差大于450m，检测概率都为1时，改进型RAIM也可以提供更好的检测性能，这也证明了改进的有效性；

对于ModeB，干扰机伪造第一号和第六号卫星的信号，同时精确控制其伪距，引导接收将第五号卫星识别为干扰；

$\left|\mathrm{COV}(({0.65s}_k^1+s_k^6)E,s_k^{5}E)\right|=0.9979---(4-48)$

因为(0.65s¹+s⁶)和之间的线性相关系数很高，接收机无法分辨y^max是被第一号卫星和第六号卫星影响还是被第五号卫星影响；图3给出了双星干扰(ModeB)的检测概率P_I和错误将第五号卫星识别为干扰的概率P_MI；从图中可以看出，传统RAIM的误识别概率P_MI随着定位误差的增大而增加，这是因为接收机首先假设只有一颗干扰卫星，而(0.65s¹+s⁶)和之间的高相似性会使接收机将第五号卫星误识别干扰星；而本文提出的改进型RAIM会对识别出的疑似干扰卫星进行验证，从而有效地降低了误识别概率，增加正确识别概率；

图4给出了双星干扰(ModeC)条件下的干扰识别概率P_I和将第六号和第七号卫星误识别为干扰的概率P_MI；因为

|COV((s⁴-0.7s⁵)E,(2.83s⁶+0.41s⁷)E)|＝1     (4-49)

因为(s⁴-0.7s⁵)和(2.83s⁶+0.41s⁷)的线性相关系数为1，那么传统RAIM完全无法区分y^max是被第四号卫星和第五号卫星影响还是被第六号卫星和第七号卫星影响，所以正确识别概率P_I和误识别概率P_MI都应该为0.5左右，图4仿真结果的也验证了这一结论；同时，改进型RAIM的正确识别概率随着定位误差的增加而增加，识别概率高于传统RAIM，这也证明了改进的有效性。

Claims (6)

1.一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，其特征在于一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法具体是按照以下步骤进行的：

步骤一、接收机状态的估计值为表示接收机三维坐标的估计值，表示接收机时钟偏差的估计值，在的邻域内得到M个粒子，组成集合其中x^(m)＝(x_U,m,y_U,m,z_U,m,δt_U,m)表示第m个接收机状态粒子，(x_U,m,y_U,m,z_U,m)表示接收机的三维坐标，δt_U,m表示接收机时钟偏差，m＝1，2....M；

步骤二、设当前可见卫星数为N，测量每一颗卫星与接收机的伪距，得到第n颗卫星与接收机的测量伪距ρ_n(n＝1,2...N)，以及卫星的三维坐标(x_n,y_n,z_n)；

步骤三、根据步骤一得到的第m个接收机状态粒子x^(m)，及步骤二得到的第n颗卫星的三维坐标(x_n,y_n,z_n)，得到第n颗卫星与第m个粒子的计算伪距

其中，c代表光速；

步骤四、根据步骤二得到的第n颗卫星与接收机的测量伪距ρ_n(n＝1,2...N)，及步骤三得到的第n颗卫星与第m个粒子的计算伪距计算第m个粒子的观测向量y^(m)＝[e_1,m,...,e_n,m,...,e_N,m]^T，其中e_n,m表示第m个粒子与第n颗卫星的伪距残差：

步骤五、根据步骤四得到的第m个粒子的观测向量y^(m)，计算第m个粒子的未归一化加权值w^(m)：

$w^{\left(m\right)}=\frac{1}{M}{\left(\frac{1}{\sqrt{2{\mathrm{\&pi;\&sigma;}}^2}}\right)}^N\exp [-\frac{\left|\right|y^{\left(m\right)}\left|\right|}{{2\mathrm{\&sigma;}}^2}];$

其中，σ表示无噪声条件下e_n,m的标准差，通常取5.9；

步骤六、根据步骤五得到的第m个粒子的未归一化加权值w^(m)，选取最大的粒子加权值w^max：

w^max＝max(w^(m))，

最大的粒子加权值w^max分别与最大观测向量y^max，最大接收机状态粒子x^max，对应的最大粒子偏差d^max相对应，其中，

$d^{\max }=\hat{x}-x^{\max },$

其中表示接收机状态的估计值；

步骤七、根据步骤六得到的最大观测向量y^max，计算附加伪距的估计值

步骤八、根据步骤七得到的附加伪距的估计值计算最大可见星联合斜率SL^max：

步骤九、根据步骤八得到的最大可见星联合斜率SL^max，计算判决门限γ；

步骤十、根据步骤九得到的判决门限γ和步骤六得到的最大观测向量y^max，若‖y^max‖≥γ则说明存在欺骗卫星，欺骗卫星标志F＝1，若不存在则跳到步骤十二，欺骗卫星标志F＝0；

步骤十一、则根据步骤六得到的最大观测向量为y^max、步骤七得到的附加伪距的估计值和残余矢量矩阵S，估计干扰卫星的序号

步骤十二、若欺骗卫星标志F＝1，去掉步骤十一的估计干扰卫星的序号的干扰，然后进行定位解算，若欺骗卫星标志F≠1直接进行定位解算；即完成了一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法。

2.根据权利要求1所述一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，其特征在于：步骤七中根据步骤六得到的最大观测向量y^max，计算附加伪距的估计值为：

$\hat{b}={\left(S^{T}S\right)}^{-1}{S}^T{y}^{\max },$

其中S＝(I-G(G^TG)^-1G^T)，为残余矢量矩阵，G代表卫星单位方向矢量矩阵，由下式得到：

$G=\left[\begin{array}{cc}{1}_1& 1\\ 1_2& 1\\ .& .\\ .& .\\ .& .\\ 1_N& 1\end{array}\right],$

1_n(n＝1,2,...,N)是地心指向第n颗卫星的单位方向矢量，I为4×4的单位矩阵。

3.根据权利要求1所述一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，其特征在于：步骤八中根据步骤七得到的附加伪距的估计值计算最大可见星联合斜率SL^max为：

${\mathrm{SL}}^{\max }=\underset{n_1,n_2,...,n_{\mathrm{Ns}}}{\max }\left[\frac{\left|\right|{\left(d^{\max }\right)}_{\mathrm{sub}}^{\mathrm{1,2,3}}\left|\right|}{\left|\right|S_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}{\hat{b}}_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}\left|\right|}\right],$

其中，1≤n₁<n₂<...<n_Ns≤N，N代表可见卫星总数，Ns表示干扰卫星总数，n₁,n₂,...,n_Ns均为整数代表可见卫星的编号，表示d^max的前三个元素构成的子矩阵，表示S的n₁,n₂,...,n_Ns列构成的子矩阵，表示的n₁,n₂,...,n_Ns行构成的子矩阵。

4.根据权利要求1所述一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，其特征在于：步骤九中根据步骤八得到的最大可见星联合斜率SL^max，计算判决门限γ为：

$\mathrm{\&gamma;}=\frac{\mathrm{PL}}{{\mathrm{SL}}^{\max }},$

其中PL为定位误差保护门限。

5.根据权利要求1所述一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，其特征在于：步骤十一中根据步骤六得到的最大观测向量为y^max、步骤七得到的附加伪距的估计值和残余矢量矩阵S，估计干扰卫星的序号为：

$({\hat{s}}_1,{\hat{s}}_2,...{\hat{s}}_{\mathrm{Ns}})=\arg \underset{n_1,n_2,...,n_{\mathrm{Ns}}}{\max }\left[\frac{\left|\right|y^{\max }\&CenterDot;\left(S_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}{\hat{b}}_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}\right)\left|\right|}{\left|\right|y^{\max }\left|\right|\left|\right|\left(S_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}{\hat{b}}_{\mathrm{sub}}^{n_1,n_2,...,n_{\mathrm{Ns}}}\right)\left|\right|}\right].$

6.根据权利要求1所述一种基于粒子滤波的改进型RAIM抗欺骗式干扰方法，其特征在于：步骤十二中定位解算具体过程为：

去掉欺骗卫星后，假设剩余的卫星数目为M，设它们的编号为：s₁,s₂......,s_M(1≤1s₁<s₂<...<s_M≤N)，接收机的坐标为(x_u,y_u,z_u)，编号为s₁的卫星与接收机之间的测量伪距为则对该卫星存在如下的伪距方程：

${\mathrm{\&rho;}}_{s_1}=\sqrt{{(x_u-x_{s_1})}^2+{(y_u-y_{s_1})}^2+{(z_u-z_{s_1})}^2}+c\&times;(\mathrm{\&Delta;T}-{\mathrm{\&Delta;t}}_{s_1})$

其中，代表编号为s₁的卫星的三维坐标，代表该卫星时钟相对于GPS时系的偏差，c代表光速；ΔT为接收机时钟相对GPS时系的偏差，令b_u＝c×ΔT，则

${\mathrm{\&rho;}}_{s_1}=\sqrt{{(x_u-x_{s_1})}^2+{(y_u-y_{s_1})}^2+{(z_u-z_{s_1})}^2}+b_u-c{\&times;\mathrm{\&Delta;t}}_{s_1}$

式中共有四个未知数，即接收机的三维坐标(x_u,y_u,z_u)和时钟偏差造成的伪距偏差b_u；对另外M-1颗卫星列伪距方程，即得到M个伪距方程，求解出接收机的三维坐标(x_u,y_u,z_u)和时钟偏差造成的伪距偏差b_u，得到接收机状态。