CN103969667B

CN103969667B - 一种基于分块粒子滤波的gnss解算方法

Info

Publication number: CN103969667B
Application number: CN201410228162.0A
Authority: CN
Inventors: 韩帅; 李亚添; 罗德巳; 孟维晓
Original assignee: Harbin Institute of Technology
Current assignee: Harbin Institute of Technology
Priority date: 2014-05-27
Filing date: 2014-05-27
Publication date: 2016-08-17
Anticipated expiration: 2034-05-27
Also published as: CN103969667A

Abstract

一种基于分块粒子滤波的GNSS解算方法，本发明涉及分块粒子滤波的GNSS解算方法。本发明是要解决传统线性化的方法计算量复杂，在迭代过程中迭代向量的维数较高，制约了在定位和定速解算的精度的问题，而提出的一种基于分块粒子滤波的GNSS解算方法。该方法是通过1、设定0时刻的和2、计算和3求归一化权值4、M_eff大于M_th，跳到5；如果M_eff<M_th，则重采样，5、求定位估计结果6、确定7、求得归一化权值8、大于跳到9；如果重采样；9、第k时刻定速解算估计结果10、k＝k+1，跳转到2，重复2到10的过程。等步骤实现的。本发明应用于分块粒子滤波的GNSS领域。

Description

一种基于分块粒子滤波的GNSS解算方法

技术领域

本发明涉及一种基于分块粒子滤波的GNSS解算方法。

背景技术

随着电子信息技术的蓬勃发展，全球导航卫星系统(Global NavigationSatellite System，GNSS)在基于位置服务中的核心作用愈发显著，其可以在全球范围内提供全天候的定位、定速和授时服务，应用领域包括飞行器导航，精密授时，大气监测，地质探测等。由于硬件技术和理论的不断进步，人们亦对GNSS系统在定位、定速解算的定位精度有了更高的要求。

接收机根据测得的伪距和伪距变化率，以及从星历中获知的卫星位置信息和速度信息，可以建立一个关于自身位置和速度信息的非线性方程组，用以定位和定速解算。最小二乘(Least Square，LS)算法是GNSS定位、定速解算的经典算法。它以最小化误差平方和为目标函数，借助牛顿迭代使非线性方程组线性化，然后利用最小二乘法实现对位置与速度的估计。基于经典最小二乘法，学者又提出了卡尔曼滤波、粒子滤波等解算算法，它们通过权值分配、平滑滤波等技术，进一步提高了定位精度。

在线性系统中，卡尔曼滤波器给出了最小均方误差意义下的最优估计。然而在实际工程运用中，大部分系统都是非线性的，这促使了非线性滤波理论的发展。比较简单的非线性滤波方法即是借助泰勒展开，将非线性系统在局部线性化，得到近似的线性系统方程，然后再利用卡尔曼滤波器实现滤波。然而线性化的方法在复杂恶劣的环境中容易发散，这限制了其应用的范围。

粒子滤波是一种基于蒙特卡洛方法的最优的估计，它通过序贯重要性采样，可以准确估计非线性非高斯场景中的状态矢量,从而实现递推贝叶斯滤波，是一种非参数化的近似最优非线性滤波方法。它不对系统方程做任何处理，同时也不需要对条件概率密度的参数做任何估计，具有较强的鲁棒性和较高的精度。因此，粒子滤波算法被应用于定位，导航和跟踪过程中。在那之后，针对基于粒子滤波的GNSS接收机结构研究逐渐增多，一些旨在提高定位精度、多径消除的改进型的方案也被提出。而现有的粒子滤波在GNSS解算的应用是单纯的是将粒子滤波用于速度和位置的整体解算，不仅计算量复杂，而且由于迭代过程中迭代向量的维数较高，制约了在定位解算和定速解算的精度。

发明内容

本发明是为了解决传统线性化的方法在复杂恶劣的环境中容易发散，计算量复杂，在迭代过程中迭代向量的维数较高，制约了在定位解算和定速解算的精度的问题，而提出了一种基于分块粒子滤波的GNSS解算方法。

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

步骤一，根据第k＝0时刻PVT解算的结果得到第k＝0时刻位置解算结果x₀和速度解算结果设定0时刻的速度矢量估计值位置矢量的估计值设定0时刻每个速度粒子权值和位置粒子权值令并且确定在0时刻M个四维位置粒子和M个四维速度粒子其中均匀分布在周围，均匀分布在的周围；定义第m个位置粒子为其中分别代表第m个位置粒子在0时刻的三维坐标，代表第m个位置粒子在0时刻状态与卫星的时钟偏差；定义第0时刻第n个卫星位置为(x_n,0,y_n,0,z_n,0)(n＝1,2,...,N)；第0时刻第m个速度粒子为：其中，分别代表第m个速度粒子在0时刻的三维坐标的速度结算结果；代表第m个速度粒子在0时刻状态与卫星的时钟偏差速度结算结果；并设置时刻的初值k＝1；转到步骤二；

步骤二、根据

\{\begin{matrix} x_{k}^{(m)} = x_{k - 1}^{(m)} + {\overset{\hat{\cdot}}{x}}_{k - 1} + e_{k} \\ y_{k}^{(m)} = h 1 (x_{k}^{(m)}, ρ_{k}) + ϵ_{k} \end{matrix}, k = 1,2 . . . . n,

计算和如果k＝1，那么和如果k>1，那么和分别来自上次迭代的步骤九和步骤二；为第k时刻第m个位置粒子的解算信息，在k时刻所有M个四维位置粒子为为第k-1时刻第m个速度粒子的定速估计结果；其中，e_k为第k时刻定位解算的过程噪声和ε_k为第k时刻定位解算的观测噪声，ρ_k为接收机得到的第k时刻的伪距值；h1(·)为定位解算观测方程，为第k时刻第m个位置粒子解算的N维观矢量；根据步骤二可以确定然后进行步骤三；

步骤三、基于过程噪声和观测噪声为不相关的高斯白噪声，将和带入公式其中为第k时刻第m个位置粒子N维观测矢量中的第n个元素；当k＝1，当k>1，上次迭代的步骤三；求出k时刻第m个位置粒子的未归一化权重根据所有第k时刻所有M个位置粒子权值和公式求得第k时刻，第m个位置粒子的归一化权值其中，σ²为过程噪声的方差；

步骤四、根据步骤三计算得出的计算第k时刻位置粒子有效粒子数量如果M_eff大于设定的有效位置粒子门限M_th，跳到步骤五；如果M_eff<M_th，则进行重采样，从中随机挑选出M_loop个粒子，组成新的粒子集合，在M_loop个粒子中，第m个位置粒子被挑中的概率为在重采样之后，每个随机挑选出的粒子的权值为1/M，其中，M_loop为重采样随机采的粒子数目；

步骤五、在对权值进行归一化之后，根据公式求得k时刻的定位估计结果

步骤六、根据公式

\{\begin{matrix} {\overset{\cdot}{x}}_{k}^{(m)} = {\overset{\cdot}{x}}_{k - 1}^{(m)} + {\overset{\cdot}{e}}_{k} \\ {\overset{\cdot}{y}}_{k}^{(m)} = h 2 ({\hat{x}}_{k}, {\overset{\cdot}{x}}_{k}^{(m)}, ρ_{k}, {\overset{\cdot}{ρ}}_{k}) + {\overset{\cdot}{ϵ}}_{k} \end{matrix}

确定如果k＝1的，那么如果k>1，来自上一次迭代的步骤六；其中，ρ_k为接收机获得的第k时刻的伪距，为接收机获得的第k时刻的伪距变化量，为定速解算的过程噪声，为定速解算的观测噪声，h2(·)为第k时刻第m个速度粒子解算的观测方程；即第k时刻第m个速度粒子的解算信息；为第k时刻第m个速度粒子的N维观测失量，为k时刻的第m个位置粒子定位结果；h2(·)定速解算的观测方程；

步骤七、基于过程噪声和观测噪声为不相关的高斯白噪声，根据和公式求第k时刻第m个速度粒子的未归一化权重其中，如果k＝1,那么如果k>1，那么来自上次迭代的步骤七；求得所有M个速度粒子未归一化权重后；根据公式归一化为第k时刻，第m个速度粒子的归一化权值其中，是的第n个元素，来自步骤六；为速度解算中的过程噪声的方差；

步骤八、根据步骤七得出的计算第k时刻有效速度粒子数量如果大于设定的有效粒子门限跳到步骤九；如果则进行重采样，从中随机挑选出个粒子，组成新的粒子集合，第m个粒子被挑中的概率为在重采样之后，每个随机挑选出的粒子的权值为1/M；然后进行步骤九，其中，为重采样随机采的速度粒子数目；

步骤九、在对权值进行归一化之后，根据和公式得到第k时刻定速解算估计结果然后进行步骤十；

步骤十，根据步骤二到步骤九求得第k时刻的定位估计结果和定速估计结果然后令k＝k+1，跳转到步骤二，重复步骤二到步骤十的过程，进行迭代；即完成了一种基于分块粒子滤波的GNSS解算方法。

发明效果

本发明基于传统的粒子滤波算法，提出一种改进型的分块粒子滤波解算方法。该方法包含两个相互关联的子粒子滤波器，分别用于定位解算和定速解算。分块使得粒子由八维矢量变为四维矢量，在粒子数一定的条件下，低维矢量可以更好的覆盖状态矢量空间，提供更准确的估计结果粒子滤波的理论基础以及基于粒子滤波的导航解算，并且提出一种改进型的分块粒子滤波解算方法。分块粒子滤波分离了位置和速度估计，降低了状态矢量和粒子维度，可以进一步提高定位和定速精度。

分块粒子滤波相比传统粒子滤波GNSS定位方法，分块粒子滤波将定位解算和定速解算进行分块。在粒子的设定中，从传统的8维矢量，分别分解成4维的速度粒子和4维的位置粒子，减少了粒子的维数。如果粒子数目一定，在定位解算和定速解算同时并行进行的条件下，将会比传统的粒子滤波具有更好的定位精度，。总体来讲，分离了位置和速度估计两个子滤波器，这两个子滤波器并不独立，而是相互辅助与迭代。

从图3和图4可以看到，分块粒子滤波解算的定位、定速误差都小于传统粒子滤波，这也证明了所提出的算法的有效性。表3比较了不同粒子数量时，两种粒子滤波解算算法的平均误差。从表中可以看出，无论粒子数为多少，分块滤波性能都优于传统粒子滤波；传统粒子滤波和分块粒子滤波的解算误差都随着粒子数的增加而降低。粒子滤波器的估算结果会逐渐趋于最优。因此，本发明对粒子滤波进行改进，提出分块粒子滤波算法，旨在提高其定位，定速精度。

附图说明

图1是具体实施方式一提出的一种基于分块粒子滤波的GNSS解算方法流程图；

图2是实施例提出的分块粒子滤波解算流程图；

图3是具体实施方式一提出的最小二乘法、传统粒子滤波和分块粒子滤波解算定位误差比较示意图；

图4是具体实施方式一提出的最小二乘法、传统粒子滤波和分块粒子滤波解算定速误差比较示意图。

具体实施方式

具体实施方式一：本实施方式的一种基于分块粒子滤波的GNSS解算方法，具体是按照以下步骤进行的：

步骤一，根据第k＝0时刻PVT(Position,Velocity,Time)解算的结果得到第k＝0时刻位置解算结果x₀和速度解算结果设定0时刻的速度矢量估计值位置矢量的估计值设定0时刻每个速度粒子权值和位置粒子权值令并且确定在0时刻M个四维位置粒子和M个四维速度粒子其中均匀分布在周围，均匀分布在的周围；定义第m个位置粒子为其中分别代表第m个位置粒子在0时刻的三维坐标，代表第m个位置粒子在0时刻状态与卫星的时钟偏差；定义第0时刻第n个卫星位置为(x_n,0,y_n,0,z_n,0)(n＝1,2,...,N)；第0时刻第m个速度粒子为：其中，分别代表第m个速度粒子在0时刻的三维坐标的速度结算结果；代表第m个速度粒子在0时刻状态与卫星的时钟偏差速度结算结果；并设置时刻的初值k＝1；转到步骤二；

步骤二、根据

\{\begin{matrix} x_{k}^{(m)} = x_{k - 1}^{(m)} + {\overset{\hat{\cdot}}{x}}_{k - 1} + e_{k} \\ y_{k}^{(m)} = h 1 (x_{k}^{(m)}, ρ_{k}) + ϵ_{k} \end{matrix}, k = 12 . . . . n,

步骤四、根据步骤三计算得出的计算第k时刻位置粒子有效粒子数量如果M_eff大于设定的有效位置粒子门限M_th(1≤M_th≤M它可以根据用户的需要自由设定)，那么不进行重采样，跳到步骤五；如果M_eff<M_th，则进行重采样，从中随机挑选出M_loop个粒子，组成新的粒子集合，在M_loop个粒子中，第m个位置粒子被挑中的概率为在重采样之后，每个随机挑选出的粒子的权值为1/M，其中，M_loop为重采样随机采的粒子数目；

步骤五、在对权值进行归一化之后，根据公式求得k时刻的定位估计结果然后进行步骤六；

步骤六、根据公式

\{\begin{matrix} {\overset{\cdot}{x}}_{k}^{(m)} = {\overset{\cdot}{x}}_{k - 1}^{(m)} + {\overset{\cdot}{e}}_{k} \\ {\overset{\cdot}{y}}_{k}^{(m)} = h 2 ({\hat{x}}_{k}, {\overset{\cdot}{x}}_{k}^{(m)}, ρ_{k}, {\overset{\cdot}{ρ}}_{k}) + {\overset{\cdot}{ϵ}}_{k} \end{matrix}

步骤八、根据步骤七得出的计算第k时刻有效速度粒子数量如果大于设定的有效粒子门限(它可以根据用户的需要自由设定)，那么不进行重采样，跳到步骤九；如果则需要进行重采样，从中随机挑选出个粒子，组成新的粒子集合，第m个粒子被挑中的概率为在重采样之后，每个随机挑选出的粒子的权值为1/M；然后进行步骤九，其中，为重采样随机采的速度粒子数目；

步骤十，根据步骤二到步骤九求得第k时刻的定位估计结果和定速估计结果然后令k＝k+1，跳转到步骤二，重复步骤二到步骤十的过程，进行迭代如图1；即完成了一种基于分块粒子滤波的GNSS解算方法。

本实施方式效果：

本实施方式基于传统的粒子滤波算法，提出一种改进型的分块粒子滤波解算方法。该方法包含两个相互关联的子粒子滤波器，分别用于定位解算和定速解算。分块使得粒子由八维矢量变为四维矢量，在粒子数一定的条件下，低维矢量可以更好的覆盖状态矢量空间，提供更准确的估计结果粒子滤波的理论基础以及基于粒子滤波的导航解算，并且提出一种改进型的分块粒子滤波解算方法。分块粒子滤波分离了位置和速度估计，降低了状态矢量和粒子维度，可以进一步提高定位和定速精度。

从图3和图4可以看到，分块粒子滤波解算的定位、定速误差都小于传统粒子滤波，这也证明了所提出的算法的有效性。表3比较了不同粒子数量时，两种粒子滤波解算算法的平均误差。从表中可以看出，无论粒子数为多少，分块滤波性能都优于传统粒子滤波；传统粒子滤波和分块粒子滤波的解算误差都随着粒子数的增加而降低。粒子滤波器的估算结果会逐渐趋于最优。因此，本实施方式对粒子滤波进行改进，提出分块粒子滤波算法，旨在提高其定位，定速精度。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤二中为第k时刻，第m个位置粒子的N维观测矢量定义为：

y_{k}^{(m)} = [\begin{matrix} ρ_{1, k} - \sqrt{{(x_{1, k} - x_{U, k}^{(m)})}^{2} + {(y_{1, k} - y_{U, k}^{(m)})}^{2} + {(z_{1, k} - z_{U, k}^{(m)})}^{2}} - {δt}_{U, k}^{(m)} \\ ρ_{2, k} - \sqrt{{(x_{2, k} - x_{U, k}^{(m)})}^{2} + {(y_{2, k} - y_{U, k}^{(m)})}^{2} + {(z_{2, k} - z_{U, k}^{(m)})}^{2}} - {δt}_{U, k}^{(m)} \\ \cdot \cdot \cdot \\ ρ_{N, k} - \sqrt{{(x_{N, k} - x_{U, k}^{(m)})}^{2} + {(y_{N, k} - y_{U, k}^{(m)})}^{2} + {(z_{N, k} - z_{U, k}^{(m)})}^{2}} - {δt}_{U, k}^{(m)} \end{matrix}]

其中，(x_n,k,y_n,k,z_n,k)(n＝1,2,...,N)代表在第n颗卫星在k个时刻的位置坐标；为在第k时刻的第m个位置状态粒子三维位置，为相对于第n颗卫星的时钟偏差。其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：ρ_k为N维伪距矢量，定义为：

ρ_k＝[ρ_1,k,ρ_2,k,...,ρ_N,k]^T (5)

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

ρ_{n, k} = \sqrt{{(x_{n, k} - x_{U, k})}^{2} + {(y_{n, k} - y_{U, k})}^{2} + {(z_{n, k} - z_{U, k})}^{2}} + c \cdot {δt}_{U, 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颗卫星的时钟偏差，c为光速。其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤五为观测方程即N维观测矢量，定义为：

其中，x_N,ky_N,kz_N,k为第n颗卫星在第k个时刻的位置坐标为第n颗卫星在第k个时刻的第m个位置状态粒子速率坐标在单位时隙中的变化量；为在第k时刻的第m个位置状态粒子速率三维位置在单位时隙中的变化量，为在第k时刻的第m个位置状态粒子三维位置的估计值，为相对于第n颗卫星的第k时刻的第m个位置状态粒子速率的时钟偏差。其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：为N维伪距变化矢量，定义为：

{\overset{\cdot}{ρ}}_{k} = {[{\overset{\cdot}{ρ}}_{1, k}, {\overset{\cdot}{ρ}}_{2, k}, . . ., {\overset{\cdot}{ρ}}_{N, k}]}^{T}

{\overset{\cdot}{ρ}}_{n, k} = \sqrt{{({\overset{\cdot}{x}}_{n, k} - {\overset{\cdot}{x}}_{U, k})}^{2} + {({\overset{\cdot}{y}}_{n, k} - {\overset{\cdot}{y}}_{U, k})}^{2} + {({\overset{\cdot}{z}}_{n, k} - {\overset{\cdot}{z}}_{U, k})}^{2}} + c \cdot {\overset{\cdot}{δ} t}_{U, k}

为第n颗卫星在第k个时刻的第m个位置状态粒子速率坐标在单位时隙中的变化量；为在第k时刻的第m个位置状态粒子速率三维位置在单位时隙中的变化量，为相对于第n颗卫星的第k时刻的第m个位置状态粒子速率的时钟偏差，c为光速。其它步骤及参数与具体实施方式一至四之一相同。

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

实施例一：

本实施例一种基于分块粒子滤波的GNSS解算方法，具体是按照以下步骤制备的：

步骤一，根据第k＝0时刻PVT(Position,Velocity,Time)解算的结果得到第k＝0时刻位置解算结果x₀＝(-2654200m，3564100m，4565300m)和速度解算结果设定0时刻的速度矢量估计值位置矢量的估计值设定0时刻每个速度粒子和位置粒子权值令并且确定在0时刻M(M＝100)个四维位置粒子和M个四维速度粒子其中均匀分布在周围，均匀分布在的周围；定义第m个位置粒子为其中分别代表第m个位置粒子在0时刻的三维坐标，代表第m个位置粒子在0时刻状态与卫星的时钟偏差；定义第0时刻第n个卫星位置为(x_n,0,y_n,0,z_n,0)(n＝1,2,...,N),N取8；假设当前接收机的可见星为8颗，可见星的位置坐标如表1所示；第0时刻第m个速度粒子为：其中，分别代表第m个速度粒子在0时刻的三维坐标的速度结算结果；代表第m个速度粒子在0时刻状态与卫星的时钟偏差速度结算结果；并设置时刻的初值k＝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

表1 可见星的位置坐标

步骤二、根据

\{\begin{matrix} x_{k}^{(m)} = x_{k - 1}^{(m)} + {\overset{\hat{\cdot}}{x}}_{k - 1} + e_{k} \\ y_{k}^{(m)} = h 1 (x_{k}^{(m)}, ρ_{k}) + ϵ_{k} \end{matrix}, k = 12 . . . . n,

步骤三、基于过程噪声和观测噪声为不相关的高斯白噪声，将和带入公式其中为第k时刻第m个位置粒子N维观测矢量中的第n个元素；当k＝1，当k>1，上次迭代的步骤三；求出k时刻第m个位置粒子的未归一化权重根据所有第k时刻所有M个位置粒子权值和公式求得第k时刻，第m个位置粒子的归一化权值其中，σ²为过程噪声的方差，设定σ＝5.9m；

步骤四、根据步骤三计算得出的计算第k时刻位置粒子有效粒子数量如果M_eff大于设定的有效位置粒子门限M_th＝2M/3(1≤M_th≤M)，那么不进行重采样，跳到步骤五；如果M_eff<M_th，则进行重采样，从中随机挑选出M_loop个粒子，组成新的粒子集合，在M_loop＝M个粒子中，第m个位置粒子被挑中的概率为在重采样之后，每个随机挑选出的粒子的权值为1/M，其中，M_loop为重采样随机采的粒子数目；

步骤六、根据公式

\{\begin{matrix} {\overset{\cdot}{x}}_{k}^{(m)} = {\overset{\cdot}{x}}_{k - 1}^{(m)} + {\overset{\cdot}{e}}_{k} \\ {\overset{\cdot}{y}}_{k}^{(m)} = h 2 ({\hat{x}}_{k}, {\overset{\cdot}{x}}_{k}^{(m)}, ρ_{k}, {\overset{\cdot}{ρ}}_{k}) + {\overset{\cdot}{ϵ}}_{k} \end{matrix}

步骤七、基于过程噪声和观测噪声为不相关的高斯白噪声，根据和公式求第k时刻第m个速度粒子的未归一化权重其中，如果k＝1,那么如果k>1，那么来自上次迭代的步骤七；求得所有M个速度粒子未归一化权重后；根据公式归一化为第k时刻，第m个速度粒子的归一化权值其中，是的第n个元素，来自步骤六；为速度解算中的过程噪声的方差,设定σ＝0.9m/s；

步骤八、根据步骤七得出的计算第k时刻有效速度粒子数量如果大于设定的有效粒子门限那么不进行重采样，跳到步骤九；如果则需要进行重采样，从中随机挑选出个粒子，组成新的粒子集合，第m个粒子被挑中的概率为在重采样之后，每个随机挑选出的粒子的权值为1/M；然后进行步骤九，其中，为重采样随机采的速度粒子数目；

步骤十，根据步骤二到步骤九求得第k时刻的定位估计结果和定速估计结果然后令k＝k+1，跳转到步骤二，重复步骤二到步骤十的过程，进行迭代；对分块粒子滤波解算算法进行仿真验证，并将其性能与传统粒子滤波和传统最小二乘法进行比较结果如下表所示：

表2最小二乘法、传统粒子滤波和分块粒子滤波解算的平均误差

从图3和图4可以看到，分块粒子滤波解算的定位、定速误差都小于传统粒子滤波，这也证明了所提出的算法的有效性；表3比较了不同粒子数量时，两种粒子滤波解算算法的平均误差；从表中可以看出，无论粒子数为多少，分块滤波性能都优于传统粒子滤波；传统粒子滤波和分块粒子滤波的解算误差都随着粒子数的增加而降低；粒子滤波器的估算结果会逐渐趋于最优。

Claims

1.一种基于分块粒子滤波的GNSS解算方法，其特征在于：一种基于分块粒子滤波的GNSS解算方法具体是按照以下步骤进行的：

步骤一，根据第k＝0时刻PVT解算的结果得到第k＝0时刻位置解算结果x₀和速度解算结果设定0时刻的速度矢量估计值位置矢量的估计值设定0时刻每个速度粒子权值和位置粒子权值令并且确定在0时刻M个四维位置粒子和M个四维速度粒子其中均匀分布在周围，均匀分布在的周围；定义第m个位置粒子为其中分别代表第m个位置粒子在0时刻的三维坐标，代表第m个位置粒子在0时刻状态与卫星的时钟偏差；定义第0时刻第n个卫星位置为(x_n,0,y_n,0,z_n,0)，n＝1,2,...,N；第0时刻第m个速度粒子为：其中，分别代表第m个速度粒子在0时刻的三维坐标的速度结算结果；代表第m个速度粒子在0时刻状态与卫星的时钟偏差速度结算结果；并设置时刻的初值k＝1；转到步骤二；

步骤二、根据计算和为第k时刻第m个位置粒子的解算信息，在k时刻所有M个四维位置粒子为为第k-1时刻第m个速度粒子的定速估计结果；其中，e_k为第k时刻定位解算的过程噪声、ε_k为第k时刻定位解算的观测噪声，ρ_k为接收机得到的第k时刻的伪距值；h1(·)为定位解算观测方程，为第k时刻第m个位置粒子解算的N维观矢量；

步骤三、基于过程噪声和观测噪声为不相关的高斯白噪声，将和带入公式其中为第k时刻第m个位置粒子N维观测矢量中的第n个元素；求出k时刻第m个位置粒子的未归一化权重根据所有第k时刻所有M个位置粒子权值和公式求得第k时刻，第m个位置粒子的归一化权值其中，σ²为过程噪声的方差；

步骤六、根据公式确定其中，ρ_k为接收机获得的第k时刻的伪距，为接收机获得的第k时刻的伪距变化量，为定速解算的过程噪声，为定速解算的观测噪声，即第k时刻第m个速度粒子的解算信息；为第k时刻第m个速度粒子的N维观测失量，h2(·)定速解算的观测方程；

步骤七、基于过程噪声和观测噪声为不相关的高斯白噪声，根据和公式求第k时刻第m个速度粒子的未归一化权重求得所有M个速度粒子未归一化权重后；根据公式归一化为第k时刻，第m个速度粒子的归一化权值其中，是的第n个元素，来自步骤六；为速度解算中的过程噪声的方差；

2.根据权利要求1所述一种基于分块粒子滤波的GNSS解算方法，其特征在于：步骤二中为第k时刻，第m个位置粒子的N维观测矢量定义为：

y_{k}^{(m)} = [\begin{matrix} ρ_{1, k} - \sqrt{{(x_{1, k} - x_{U, k}^{(m)})}^{2} + {(y_{1, k} - y_{U, k}^{(m)})}^{2} + {(z_{1, k} - z_{U, k}^{(m)})}^{2}} - δ t_{U, k}^{(m)} \\ ρ_{2, k} - \sqrt{{(x_{2, k} - x_{U, k}^{(m)})}^{2} + {(y_{2, k} - y_{U, k}^{(m)})}^{2} + {(z_{2, k} - z_{U, k}^{(m)})}^{2}} - {δt}_{U, k}^{(m)} \\ ... \\ ρ_{N, k} - \sqrt{{(x_{N, k} - x_{U, k}^{(m)})}^{2} + {(y_{N, k} - y_{U, k}^{(m)})}^{2} + {(z_{N, k} - z_{U, k}^{(m)})}^{2}} - δ t_{U, k}^{(m)} \end{matrix}]

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

3.根据权利要求2所述一种基于分块粒子滤波的GNSS解算方法，其特征在于：

ρ_k为N维伪距矢量，定义为:

ρ_k＝[ρ_1,k,ρ_2,k,...,ρ_N,k]^T (5)

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

ρ_{n, k} = \sqrt{{(x_{n, k} - x_{U, k})}^{2} + {(y_{n, k} - y_{U, k})}^{2} + {(z_{n, k} - z_{U, k})}^{2}} + c \cdot {δt}_{U, k}

其中，(x_U,k,y_U,k,z_U,k)为在第k时刻的状态粒子三维位置，δt_U,k为相对于第n颗卫星的时钟偏差，c为光速。

4.根据权利要求3所述一种基于分块粒子滤波的GNSS解算方法，其特征在于：步骤六为观测方程即N维观测矢量，定义为：

{\overset{\cdot}{y}}_{k}^{(m)} = [\begin{matrix} {\overset{\cdot}{ρ}}_{1, k} - [\begin{matrix} (x_{1, k} - {\hat{x}}_{U, k}) ({\overset{\cdot}{x}}_{1, k} - {\overset{\cdot}{x}}_{U, k}^{(m)}) \\ + (y_{1, k} - {\hat{y}}_{U, k}) ({\overset{\cdot}{y}}_{1, k} - {\overset{\cdot}{y}}_{U, k}^{(m)}) \\ + (z_{1, k} - {\hat{z}}_{U, k}) ({\overset{\cdot}{z}}_{1, k} - {\overset{\cdot}{z}}_{U, k}^{(m)}) \end{matrix}] / ρ_{1, k} - \overset{\cdot}{δ} t_{U, k}^{(m)} \\ {\overset{\cdot}{ρ}}_{2, k} - [\begin{matrix} (x_{2, k} - {\hat{x}}_{U, k}) ({\overset{\cdot}{x}}_{2, k} - {\overset{\cdot}{x}}_{U, k}^{(m)}) \\ + (y_{2, k} - {\hat{y}}_{U, k}) ({\overset{\cdot}{y}}_{2, k} - {\overset{\cdot}{y}}_{U, k}^{(m)}) \\ + (z_{2, k} - {\hat{z}}_{U, k}) ({\overset{\cdot}{z}}_{2, k} - {\overset{\cdot}{z}}_{U, k}^{(m)}) \end{matrix}] / ρ_{2, k} - \overset{\cdot}{δ} t_{U, k}^{(m)} \\ ... \\ {\overset{\cdot}{ρ}}_{N, k} - [\begin{matrix} (x_{N, k} - {\hat{x}}_{U, k}) ({\overset{\cdot}{x}}_{N, k} - {\overset{\cdot}{x}}_{U, k}^{(m)}) \\ + (y_{N, k} - {\hat{y}}_{U, k}) ({\overset{\cdot}{y}}_{N, k} - {\overset{\cdot}{y}}_{U, k}^{(m)}) \\ + (z_{N, k} - {\hat{z}}_{U, k}) ({\overset{\cdot}{z}}_{N, k} - {\overset{\cdot}{z}}_{U, k}^{(m)}) \end{matrix}] / ρ_{n, k} - \overset{\cdot}{δ} t_{U, k}^{(m)} \end{matrix}]

其中，为第n颗卫星在第k个时刻的第m个位置状态粒子速率坐标在单位时隙中的变化量；为在第k时刻的第m个位置状态粒子速率三维位置在单位时隙中的变化量，为在第k时刻的第m个位置状态粒子三维位置的估计值，为相对于第n颗卫星的第k时刻的第m个位置状态粒子速率的时钟偏差。

5.根据权利要求4所述一种基于分块粒子滤波的GNSS解算方法，其特征在于：为N维伪距变化矢量，定义为：

{\overset{\cdot}{ρ}}_{k} = {[{\overset{\cdot}{ρ}}_{1, k}, {\overset{\cdot}{ρ}}_{2, k}, ..., {\overset{\cdot}{ρ}}_{N, k}]}^{T}

{\overset{\cdot}{ρ}}_{n, k} = \sqrt{{({\overset{\cdot}{x}}_{n, k} - {\overset{\cdot}{x}}_{U, k})}^{2} + {({\overset{\cdot}{y}}_{n, k} - {\overset{\cdot}{y}}_{U, k})}^{2} + {({\overset{\cdot}{z}}_{n, k} - {\overset{\cdot}{z}}_{U, k})}^{2}} + c \cdot \overset{\cdot}{δ} t_{U, k} .