Particle filter resampling-based GNSS multipath correction method
-
摘要:
多路径是制约高精度卫星定位服务的重要因素之一,当前多路径模型中未能有效顾及噪声分布与处理时效性等问题. 本文提出了一种基于粒子滤波重采样的GNSS伪距多路径修正方法. 首先,利用原始伪距多路径序列构建粒子滤波初始化阶段状态空间模型;其次,顾及粒子更新过程噪声传递并预测下一时刻粒子群及其权重分布;最后,基于顾及卫星高度角的重采样技术更新多路径粒子,加权并滑动输出单历元多路径延迟量. 通过MGEX静态观测数据分析表明,本文提出的GNSS多路径修正方法分别可实现北斗卫星导航系统(BeiDou Navigation Satellite System,BDS)、GPS、GLONASS与Galileo伪距多路径标准差(standard deviation,STD)降低59.2%、58%、61.7%与59.8%;精密单点定位(precise point positioning,PPP)精度在E、N与U方向分别提升了10.68%、11.95%与26.13%;动态定位精度在E、N与U方向分别提升了5.9%、1.1%与53.4%. 因此,本文提出的修正方法在处理GNSS多路径对位置服务性能提升具有一定意义,并且对完善数据预处理模型具有重要作用.
Abstract:Multipath is a critical factor hindering the provision of high-precision satellite positioning services. Notably, current multipath models have failed to effectively address the challenges posed by noise distribution and processing timeliness. In this paper, we introduce a novel correction method for GNSS pseudorange multipath, leveraging particle filter resampling. Initially, we establish the initial state space model of the particle filter, utilizing the original pseudorange multipath sequence. Subsequently, we consider the noise propagation during the particle update process and predict the distribution of particles and their weights at the subsequent time step. Crucially, we update the multipath particles, taking into account the satellite elevation angle, to obtain a weighted and sliding output of the single epoch multipath delay. Analysis of MGEX static observation data reveals that our proposed GNSS multipath correction method significantly improves the STD by 59.2%, 58%, 61.7%, and 59.8% for BDS, GPS, GLONASS, and Galileo pseudorange multipath, respectively. The PPP positioning accuracy has undergone significant enhancements, recording improvements of 10.68%, 11.95%, and 26.13% in the E, N, and U directions, respectively. The accuracy of dynamic positioning in the E, N, and U directions underwent respective enhancements of 5.9%, 1.1%, and a substantial 53.4%. Consequently, the proposed correction method for addressing GNSS multipath effects holds substantial significance in augmenting the performance of location services and plays a pivotal role in enhancing data preprocessing models.
-
Keywords:
- GNSS /
- multipath effects /
- particle filter /
- resampling /
- observation noise
-
0. 引 言
多路径效应是GNSS误差的主要来源,无法通过差分技术有效地消除[1-4]. 此外,由于测站点周围环境的多样性,载波相位观测误差可能达到cm级,而伪距观测误差则可能达到m级[5-7]. 因此,深入研究和分析多路径误差的变化特性对于提升GNSS定位精度至关重要[8-11].
为实现GNSS高精度定位,各学者提出了诸多多路径缓解策略[12-13]. 最直接的方法就是选择低多路径环境作为测站点,但在实际应用中,这种方法缺少可行性[14-16]. 从优化接收机硬件层面考虑[17-18],使用大口径抛物面天线以增强GNSS信号的追踪能力[19]或利用天线运动来缓解多路径效应[20],确实能够取得一定的效果. 然而,这些方法受限于尺寸和安全性等问题,并不适用于移动载体. 通常情况下,观测数据都不可避免地会受到多路径效应的影响,那么就需要在数据处理过程中借助数学模型来减轻其影响. 例如,根据信噪比(signal to noise ratio,SNR)[21-22]或载噪比(carrier to noise ratio,CNR)[23]的大小对GNSS观测数据进行加权处理,以此减小多路径效应产生的误差. 尽管优化硬件和数学模型都能在一定程度上缓解多路径问题,但多路径延迟并未得到真正的校准. 甚至一些原本具有较高SNR或低高度角的观测值,即使它们并未受到多路径的影响,在经过算法处理后也可能会被降低权重[2]. 基于多路径效应的时空不变性,可以利用GNSS数据后处理中的模型进行静态环境下的实时多路径缓解;例如,恒星日滤波方法(sidereal filtering,SF)就是利用卫星运行轨道的周日重复性对多路径误差进行分离[24];在空间域方面,文献[25]根据观测值建立了基于接收机天线的多路径环境图表;文献[17]和文献[26]进一步提出了利用半天球模型(multipath hemispherical map,MHM)缓解多路径的方法,但是该方法可能受高频多路径或传输时信号衰减的限制. 此外,滤波处理技术也是缓解多路径效应的重要方法;包括Vondrak滤波[27-28]、小波分析(wavelet analysis,WA)[29]和卡尔曼滤波(Kalman filter,KF)[30];然而,这些滤波方法在处理多路径问题时,往往未能充分考虑到观测噪声对定位精度的影响. 因此,未来的研究需要进一步探索如何在减轻多路径效应的同时,更好地处理观测噪声,以实现GNSS的高精度定位.
针对GNSS伪距多路径误差与观测噪声之间的高耦合性,特别是在观测序列呈现非高斯非线性分布的情况下,粒子滤波(particle filter,PF)作为一种参数估计方法,能够很好地解决这一问题. 传统的滤波方法在处理多路径问题时,通常假设观测噪声为不相关的高斯白噪声,这使得其精度难以保证. 而粒子滤波的优势在于其不受系统模型线性或高斯假设的约束,利用样本形式描述状态概率密度,从而避免了对状态变量的概率分布做过多约束. 基于粒子滤波的优势,文献[31]使用Laplace近似得到的重要分布,有效减少了滤波所需的粒子数,降低了算法的运算量;但是随着时间的推移,该方法可能会导致粒子退化和粒子多样性的降低. 针对这一问题,文献[32]提出了一种均值漂移与粒子滤波融合算法,该算法在减少计算量的同时,有效地缓解了粒子权值的退化现象,提高了状态估计的精度. 此外,文献[33]和文献[34]在粒子滤波的框架内,通过优化粒子权重的策略,进一步提高了粒子滤波的性能. 这些优化策略不仅增强了粒子滤波的鲁棒性,还使其在高多路径环境下能够更准确地估计多路径延迟量,为GNSS伪距多路径误差的处理提供了更为有效的解决方案.
综上,本文将采用粒子滤波算法,通过运用重采样技术对粒子进行递归更新,对伪距多路径延迟量进行修正,从而提高GNSS定位精度.
1. 粒子滤波重采样方法
以北斗卫星导航系统(BeiDou Navigation Satellite System,BDS)为例,将接收机r观测到卫星s的伪距和载波相位观测方程表示如下:
$$ \begin{array}{l}\left\{\begin{array}{l}{P}_{{\mathrm{r}},i}^{{\mathrm{s}}}={\rho }_{{\mathrm{r}}}^{{\mathrm{s}}}+c\left({{\mathrm{d}}t}_{{\mathrm{r}}}+{\delta t}_{{\mathrm{r}},i}\right)-c\left({{\mathrm{d}}t}^{{\mathrm{s}}}-{\delta t}_{i}^{{\mathrm{s}}}\right)+{T}_{{\mathrm{r}}}^{{\mathrm{s}}}\\\qquad +{I}_{{\mathrm{r}},i}^{{\mathrm{s}}}+E+{M}_{{\mathrm{r}},pi}^{{\mathrm{s}}}+{\varepsilon }_{{\mathrm{r}},pi}^{{\mathrm{s}}}\\ {\lambda }_{i}{\varphi }_{{\mathrm{r}},i}^{{\mathrm{s}}}={\rho }_{{\mathrm{r}}}^{{\mathrm{s}}}+c\left({{\mathrm{d}}t}_{{\mathrm{r}}}+{\delta t}_{{\mathrm{r}},i}\right)-c\left({{\mathrm{d}}t}^{{\mathrm{s}}}-{\delta t}_{i}^{{\mathrm{s}}}\right)+{T}_{{\mathrm{r}}}^{{\mathrm{s}}}\\\qquad \;\;\;-{I}_{{\mathrm{r}},i}^{{\mathrm{s}}}+ E-{\lambda }_{i}{N}_{{\mathrm{r}},i}^{{\mathrm{s}}}+{M}_{{\mathrm{r}},\varphi i}^{{\mathrm{s}}}+{\varepsilon }_{{\mathrm{r}},\varphi i}^{{\mathrm{s}}}\end{array}\right.\end{array} $$ (1) 式中:上标s为卫星编号;下标r为接收机编号,
$ {i}({i}={1,2}) $ 表示$ \mathrm{B}1\mathrm{I} $ 和$ \mathrm{B}2\mathrm{I} $ 频率的编号;$ P $ 为伪距观测量;$ \lambda $ 为该频率的波长;$ \varphi $ 为以周为单位的载波相位观测量;$ \rho $ 为卫星s到接收机r的几何距离;$ c $ 为真空中的光速;$ {\mathrm{d}}t $ 和$ \delta t $ 分别表示钟差和硬件延迟;$ T $ 和$ I $ 分别为对流层延迟和电离层延迟;$ E $ 为卫星星历误差;$ N $ 为载波的整周模糊度;$ M $ 和$ \varepsilon $ 为表示多路径延迟量和观测噪声.为消除钟差、星历误差、大气层延迟等影响,对相同频率的伪距和载波相位观测量进行线性组合,使多路径延迟量和观测噪声准确分离;同时将不同频率下的载波相位观测量做差,获得电离层延迟
$ I $ ,即可得到波段$ \mathrm{B}1\mathrm{I} $ 和$ \mathrm{B}2\mathrm{I} $ 上的伪距多路径延迟量$$ \begin{array}{l}\left\{\begin{array}{l}{M}_{{\mathrm{r}},p1}^{{\mathrm{s}}}={P}_{{\mathrm{r}},1}^{{\mathrm{s}}}-\dfrac{{f}_{1}^{2}+{f}_{2}^{2}}{{f}_{1}^{2}-{f}_{2}^{2}}{\lambda }_{1}{\varphi }_{{\mathrm{r}},1}^{{\mathrm{s}}}+\dfrac{2{f}_{2}^{2}}{{f}_{1}^{2}-{f}_{2}^{2}}{\lambda }_{2}{\varphi }_{{\mathrm{r}},2}^{{\mathrm{s}}}-{\varepsilon }_{{\mathrm{r}},p1}^{{\mathrm{s}}}+K\\ {M}_{{\mathrm{r}},p2}^{{\mathrm{s}}}={P}_{{\mathrm{r}},2}^{{\mathrm{s}}}-\dfrac{{2f}_{1}^{2}}{{f}_{1}^{2}-{f}_{2}^{2}}{\lambda }_{1}{\varphi }_{{\mathrm{r}},1}^{{\mathrm{s}}}+\dfrac{{f}_{1}^{2}+{f}_{2}^{2}}{{f}_{1}^{2}-{f}_{2}^{2}}{\lambda }_{2}{\varphi }_{{\mathrm{r}},2}^{{\mathrm{s}}}-{\varepsilon }_{{\mathrm{r}},p2}^{{\mathrm{s}}}+K\end{array}\right.\end{array} $$ (2) 式中:
$ {f}_{1} $ 和$ {f}_{2} $ 分别为$ \mathrm{B}1\mathrm{I} $ 和$ \mathrm{B}2\mathrm{I} $ 的频率;$ K $ 为整周模糊度组合,一般为常数,不会对多路径延迟量的变化趋势产生影响. 由于载波相位观测量的精度远大于伪距观测量的精度,为了进一步简化观测模型,便于分析,已忽略式(1)中载波相位多路径延迟量和观测噪声的影响;因此,伪距观测噪声是影响伪距多路径误差大小的重要因素.基于卫星与接收机的空间关系,利用常加速(constant acceleration,CA)模型,建立伪距多路径延迟量的状态转移模型,那么系统在历元间的状态转移模型和观测模型可表示为
$$ \left\{\begin{array}{l}\begin{bmatrix}{M}_{k}\\ {M}'_{k}\\ {M}''_{k}\end{bmatrix}=\begin{bmatrix}1& \Delta t& \dfrac{\Delta {t}^{2}}{2}\\ 0& 1& \Delta t\\ 0& 0& 1\end{bmatrix}\begin{bmatrix}{M}_{k-1}\\ {M}'_{k-1}\\ {M}''_{k-1}\end{bmatrix}+\begin{bmatrix}\dfrac{\Delta {t}^{3}}{6}\\ \dfrac{\Delta {t}^{2}}{2}\\ \Delta t\end{bmatrix}{{\boldsymbol{\varPhi}} }_{k-1}\\ {{\boldsymbol{M}}}_{k}=\begin{bmatrix}1& 0& 0\end{bmatrix}\begin{bmatrix}{M}_{k}\\ {M}'_{k}\\ {M}''_{k}\end{bmatrix}+{{\boldsymbol{\psi}} }_{k}\end{array}\right. $$ (3) 式中:
$ {{\boldsymbol{M}}}_{k} $ 为历元$ {t}_{k} $ 时刻的伪距多路径延迟量,$ {M}'_{k} $ 和$ {M}''_{k} $ 分别为延迟量的变化率和加速率;$ \Delta t $ 为相邻历元时间间隔;$ {{\boldsymbol{\varPhi}} }_{k-1} $ 为系统在$ {t}_{k-1} $ 时刻产生的过程噪声;$ {{\boldsymbol{\psi}} }_{k} $ 为系统在$ {t}_{k} $ 时刻的观测噪声. 上式可表示为矩阵形式:$$ \begin{array}{l}\left\{\begin{array}{l}{\boldsymbol{X}}_{k}={\boldsymbol{FX}}_{k-1}+{\boldsymbol{B}}{{\boldsymbol{\varPhi}} }_{k-1}\\ {\boldsymbol{L}}_{k}={\boldsymbol{HX}}_{k}+{{\boldsymbol{\psi}} }_{k}\end{array}\right.\end{array} $$ (4) 式中,在历元
$ {t}_{k} $ 时刻,观测系统状态向量$ {\boldsymbol{X}}_{k}= {\left[{M}_{k},{M}'_{k},{M}''_{k}\right]}^{{\mathrm{T}}} $ ,观测向量$ {\boldsymbol{L}}_{k}={{\boldsymbol{M}}}_{k} $ ,$ \boldsymbol{F} $ 、$ \boldsymbol{H} $ 、B 分别为$ {t}_{k-1} $ 到$ {t}_{k} $ 时刻伪距多路径延迟量的状态转移矩阵、观测转移矩阵和过程噪声矩阵.设定粒子数为N,对滤波进行初始化. 开始生成N个在区间
$ \left[v\left(0\right)-{3\sigma }_{v},v\left(0\right)+{3\sigma }_{v}\right] $ 均匀分布的粒子$ {\left\{{x}_{0}^{i}\right\}}_{i=1}^{N} $ ;其中$ {\sigma }_{v} $ 为误差序列$v $ 的均方根误差(root mean squared error,RMSE),利用RMSE计算公式可以求得;同时令所有粒子的初始权重均为1/N.针对伪距多路径延迟量的非高斯非线性问题,在贝叶斯理论的预测更新迭代过程中,难以得到后验概率的解析解,因此在粒子滤波中利用蒙特卡洛采样代替计算后验概率. 利用采样粒子估计分布函数的期望值表示为
$$ E\left[f\left({x}_{n}\right)\right]\approx \frac{1}{N}\sum _{i=1}^{N}f\left({x}_{n}^{\left(i\right)}\right) $$ (5) 其中,
$ {x}_{n}^{\left(i\right)} $ 表示第$ i $ 个采样粒子;$ f\left({x}_{n}\right) $ 表示采样粒子服从的分布函数,根据重要性采样有$$E\left[f\left({x}_{k}\right)\right]=\frac{{E}_{q\left({x}_{k}|{y}_{1:k}\right)}\left[{\omega }_{k}\left({x}_{k}\right)*f\left({x}_{k}\right)\right]}{{E}_{q\left({x}_{k}|{y}_{1:k}\right)}\left[{\omega }_{k}\left({x}_{k}\right)\right]}$$ (6) 其中,
$ q\left({x}_{k}|{y}_{1:k}\right) $ 表示重要性概率密度函数;$ {\omega }_{k}\left({x}_{k}\right) $ 表示采样粒子所占权重.联合式(5)~(6)可得
$$E\left[f\left({x}_{k}\right)\right]\approx \sum _{i=1}^{N}{\widetilde{\omega }}_{k}\left({x}_{k}^{\left(i\right)}\right)*f\left({x}_{k}^{\left(i\right)}\right) $$ (7) 式(7)中,
$ {\widetilde{\omega }}_{k}\left({x}_{k}^{\left(i\right)}\right)=\frac{\displaystyle\omega_{k}\left({x}_{k}^{\left(i\right)}\right)}{\displaystyle\sum\limits _{i=1}^{N}{\omega }_{k}\left({x}_{k}^{\left(i\right)}\right)} $ ,表示归一化后的权重.根据序贯重要性采样对粒子权重进行递归
$${\omega }_{k}^{\left(i\right)}\propto {\omega }_{k-1}^{\left(i\right)}\frac{p\left({y}_{k}|{x}_{k}^{\left(i\right)}\right)p\left({x}_{k}^{\left(i\right)}|{x}_{k-1}^{\left(i\right)}\right)}{q\left({x}_{k}^{\left(i\right)}|{x}_{0:k-1}^{\left(i\right)},{y}_{1:k}\right)}$$ (8) 其中,
$ p\left({y}_{k}|{x}_{k}^{\left(i\right)}\right) $ 表示观测噪声分布,由观测方程决定,它的概率分布形状和观测值的一模一样;$ p\left({x}_{k}^{\left(i\right)}|{x}_{k-1}^{\left(i\right)}\right) $ 表示过程噪声分布,由状态转移方程决定.根据序贯重要性采样,假设
$ q\left({x}_{k}|{x}_{0:k},{y}_{1:k}\right)= q\left({x}_{k}|{x}_{k-1},{y}_{k}\right) $ ,并结合马尔科夫模型,那么粒子权重的递归形式可表示为$${\omega }_{k}^{\left(i\right)}\propto {\omega }_{k-1}^{\left(i\right)}p\left({y}_{k}|{x}_{k}^{\left(i\right)}\right) $$ (9) 由此,k时刻伪距多路径误差的最小均方差估计可表示为
$$ {\widehat{x}}_{k}\approx \sum _{i=1}^{N}\widetilde{\omega}\left({x}_{k}^{i}\right){x}_{k}^{i} $$ (10) 在采用标准的粒子滤波算法时,通常选择先验概率密度函数作为重要性密度函数,以替代后验概率密度函数来执行采样过程. 然而,这种方法存在一个明显的不足:由于先验概率密度函数在构建时并未纳入当前的伪距观测值,这可能导致从重要性密度函数中抽取的样本与真实后验概率密度函数下的样本存在显著偏差. 这种偏差会随着迭代过程的进行而逐渐累积,使得粒子的权重越来越集中于少数几个粒子,而其他粒子的权重则逐渐减小. 这种权重分布的不均匀性会严重限制新生成粒子对实际观测值后验概率分布的准确反映能力.
为了解决这一问题,可以引入重采样技术. 有效样本数
$ {N}_{e} $ 是评估粒子权重分布均匀性的关键指标,它近似表示了当前粒子集合中能有效反映后验概率分布的粒子数量. 通过监测$ {N}_{e} $ 的值,我们可以及时判断是否需要进行重采样,以保持粒子权重的均衡,从而确保粒子滤波算法能够更准确地描述系统的状态.$${\widehat{N}}_{e}\approx \frac{1}{\displaystyle\sum\limits _{i=1}^{N}{\left[\omega \left({x}_{k}^{i}\right)\right]}^{2}} $$ (11) 当有效样本数较少时,表明粒子的退化程度较为严重. 通过增加实际粒子数量可减轻粒子的退化程度. 当有效样本数
$ {N}_{e} $ 小于阈值$ {N}_{t} $ 时执行重采样操作.在进行采样前,由于卫星高度角越低时多路径效应越明显,导致误差波动量级更大,进而影响滤波效果. 因此采取了一种基于卫星高度角变化的策略来控制重采样粒子数. 首先基于高度角正弦值,根据k时刻归一化后的粒子权重集合,构建了一个多项式分布. 随后,从这个分布中抽取N个粒子权重,并相应生成N个序列号,每个序列号均为1到N之间的整数;接着定义序号值
$ {N}_{i} $ 等于$ i $ 的粒子总数,并统计每个序号值相等的粒子数量,从而得到$ {\left\{{N}_{i}\right\}}_{i=1}^{N} $ . 基于统计结果采取复制策略,将$ {\left\{{x}_{k}^{i}\right\}}_{i=1}^{N} $ 中每个粒子$ {x}_{k}^{i} $ 复制$ {N}_{i} $ 次,最终得到重采样后的粒子集合$ {\left\{{x}_{k}^{i*}\right\}}_{i=1}^{N} $ ,且重采样后的粒子集合中各粒子的权值等于1/N,确保了粒子权重的均匀分布,从而提高了滤波的准确性和鲁棒性.由此,经重采样后的k时刻伪距多路径误差的最小均方差估计可表示为
$$ {\widehat{x}}_{k}\approx \frac{1}{N}\displaystyle\sum\limits _{i=1}^{N}{x}_{k}^{i*} $$ (12) 由此可根据式(9)~(12)实现基于重采样粒子滤波的递归更新,达到伪距观测系统状态的实时估计. 具体实现流程如图1所示.
2. 伪距多路径修正实验分析
为分析本文提出的顾及粒子滤波重采样过程的GNSS多路径修正方法的效果,分别对BDS、GPS、GLONASS、Galileo的观测数据进行分析. 首先,在MGEX(Multi-GNSS Experiment)测站中选取2023年年积日第64~70天连续1周的GNSS观测数据,测站分布如图2所示.
然后分别对四系统的伪距观测方程和载波相位观测方程进行线性组合,提取出各系统选定站点的伪距多路径误差和观测噪声,利用重采样粒子滤波方法对数据进行去噪处理,估计站点的多路径误差序列. 图3以各系统原始和改正后的多路径误差序列作为对比方案;表1分别统计了实验中不同方案下多路径误差的RMSE和平均值(AVE). 结果表明,经过重采样粒子滤波处理后,测站点接收的各系统卫星数据中包含的伪距多路径误差波动量级更小,且均值更加趋向于0,符合正态分布的特性,验证了该方法所提取的伪距多路径误差符合实际结果,较好地实现了观测噪声的去除,有利于提高卫星定位的精度.
表 1 各系统伪距多路径误差统计m 卫星系统 RMSE AVE 原始 改正 原始 改正 BDS 0.361 5 0.147 6 0.010 4 0.005 1 GPS 0.609 6 0.256 0 −0.262 0 0.195 0 GLONASS 0.430 9 0.165 1 0.004 8 0.001 9 Galileo 0.526 2 0.159 3 −0.010 8 0.010 1 为了验证该方法在剔除多路径序列中观测噪声的有效性,对原始数据和经过改正的数据进行了处理,进而生成了伪距多路径残差序列. 如图4所示,经过该方法处理后的残差序列大致呈现出高斯分布的特征,这与观测噪声的分布情况相吻合. 因此,可以推断出,被剔除的残差序列可近似视为观测噪声,从而验证了该方法在剔除观测噪声方面的有效性.
经过对不同粒子总数(N = 50、200、400、600、800、
1000 )下的伪距多路径序列进行重采样粒子滤波处理,我们深入分析了该方法在多个系统下不同卫星数据中的处理效果. 具体对比了各系统下有效粒子数($ {N}_{e} $ )、RMSE和单历元多路径延迟量(M)解算时间的变化趋势. 结果表明,随着粒子总数(N)的增加,各系统的有效粒子数($ {N}_{e} $ )也呈现出正相关的增长趋势.由图5~7分析表明,随着粒子总数(N)的增加,各系统的有效粒子数(
$ {N}_{e} $ )也呈现出正相关的增长趋势. 这意味着当粒子总数增多时,系统能够更有效地模拟和逼近真实的概率分布,从而提高滤波的准确性和稳定性. 同时,我们还观察到,随着粒子总数的增加,误差序列的波动也逐渐减小. 这是因为粒子总数的增加使得滤波算法在估计过程中拥有更多的样本,能够更全面地覆盖真实状态的可能取值范围,进而减少由于样本不足导致的估计偏差和误差波动. 另一方面,尽管有效粒子数($ {N}_{e} $ )、RMSE随着粒子总数的增加而表现得更好,但处理多路径延迟量(M)时所需的时间也随之增加. 这表明,在提升滤波效果的同时,也需要权衡计算时间和计算资源的使用.为进一步验证该方法的准确性和可靠性,以6个实测点(Hzjy、Deqi、Lina、Fuya、Xish、Keqi)2018年年积日第251~255天连续5 d的静态实验观测数据为基础,通过重采样粒子滤波算法对伪距多路径误差和观测噪声进行处理. 图8表示的是各站点对误差序列滤波前后的比较,虽然各个站点周围环境不同导致多路径误差随着环境而改变,但是各站点观测数据经过处理后的多路径误差序列相较于原始序列的波动和峰值都要更小.
经过对多路径误差序列的进一步处理,我们计算了各实测点在滤波前后的RMSE和AVE. 如表2所示,经过滤波后的误差序列展现出了更为平稳的变化幅度,并且其均值更加趋近于0. 这一显著的变化说明了该方法在去除噪声方面取得了良好的效果,能够有效地减少多路径效应对定位精度的影响,从而提高了定位数据的可靠性和准确性.
表 2 各实测点伪距多路径误差统计m 实测点 RMSE AVE 原始数据 改正数据 原始数据 改正数据 Hzjy 0.334 1 0.231 0 0.208 3 −0.151 7 Deqi 0.276 5 0.181 6 −0.108 3 0.088 2 Lina 0.291 3 0.183 5 −0.085 0 0.072 4 Fuya 0.273 2 0.191 5 0.115 8 −0.106 9 Xish 0.239 7 0.168 0 0.016 5 0.007 0 Keqi 0.500 2 0.178 9 0.002 1 −0.001 4 由于上述针对多路径误差序列处理的实验中无法说明重采样粒子滤波算法对卫星定位的精度起到优化作用,进一步利用上述各站点的观测数据进行精密单点定位(precise point positioning,PPP)精度分析. 实验中以移动站与动态站间组成的差分模型为基础,基于 Inertial Explorer 8.80 (IE) 软件计算出的位置参考序列为参考[9].为具体说明本文提出的伪距多路径修正效果,分别设定了2种定位实验策略,即传统多路径处理策略(白噪声(white noise,WN)项)和粒子滤波建模修正策略. 图9中以VACS与SGOC为例,分别针对GPS单系统与GPS+BDS组合系统,利用传统多路径处理模型与粒子滤波方案,对比相应的测站E、N与U方向的定位残差. 同时,利用选择的MGEX测站1周内的定位结果取平均值,如表3所示. 定位结果表明,采用本文所提出的粒子滤波重采样多路径处理方法后,测站点E、N和U方向的定位精度分别提升了10.68%、11.95%和26.13%. 同时,该方法也基本实现了收敛时间的缩短,其中N方向的收敛速度最快. 此外,通过对比GPS单系统与GPS+BDS组合系统,定位精度和收敛时间均得到了进一步的提升,说明了多系统观测数据的引入对定位模型的性能具有促进作用.
表 3 基于不同多路径处理模型的定位精度统计方向 GPS GPS+BDS WN 粒子滤波 提升率/% WN 粒子滤波 提升率/% E 2.88 2.88 0.00 2.34 2.09 10.68 N 2.26 1.99 11.95 1.05 1.05 0.00 U 26.18 19.34 26.13 20.18 16.33 19.08 本文进一步利用实测动态数据进行验证,实验中数据采集轨迹及GNSS接收机配置如图10所示. 为了提供动态参考位置,数据采集中配备了惯性测量设备,并利用IE软件进行相对定位解算. 数据采集时间为2023年8月10日上午,采样频率为1 Hz. 为验证本文提出的基于粒子滤波重采样的GNSS多路径修正方法在实时动态定位中的效果,以IE解算的动态轨迹为参考,根据初始阶段伪距多路径延迟量的概率密度函数,分别对原始观测数据中的多路径采用WN与粒子滤波处理策略,进行动态定位实验分析.
实验中以多系统观测数据为基础,动态定位残差序列如图11所示,相应的定位效果统计如表4所示. 通过动态实验发现,由于受显著的多路径误差影响,本文提出的粒子滤波在实时数据处理中有效的分离多路径延迟量,并进行抑制处理,尤其是U方向得到了明显改善. 相较于传统WN模型而言,本文提出的多路径粒子滤波处理算法可实现动态定位E、N与U方向分别提升5.9%,1.1%与53.4%的提升效果.
表 4 基于不同多路径处理模型的动态定位残差统计多路径策略 E N U WN/m 0.051 0.020 0.641 粒子滤波/m 0.048 0.018 0.299 提升率/% 5.9 1.1 53.4 3. 结 论
本文针对当前多路径模型中存在的噪声分布处理不足和时效性挑战,提出了一种基于粒子滤波重采样的GNSS伪距多路径修正方法. 经过原始和改正后的多路径序列对比,以及对多路径残差分布、有效粒子数(
$ {N}_{e} $ )、RMSE和单历元多路径延迟量(M)解算时间的分析,得出以下结论:1) 经过该方法处理后的残差序列大致呈现出高斯分布的特征,这与观测噪声的分布情况相吻合. 因此,我们可以推断出,被剔除的残差序列可近似视为观测噪声,从而验证了该方法在剔除观测噪声方面的有效性.
2) 随着粒子总数(N)的增加,各系统的有效粒子数(
$ {N}_{e} $ )也呈现出正相关的增长趋势. 同时,误差序列的波动也逐渐减小. 这意味着当粒子总数增多时,系统能够更有效地模拟和逼近真实的概率分布,并使得滤波算法在估计过程中拥有更多的样本,能够更全面地覆盖真实状态的可能取值范围.3) 另一方面,尽管有效粒子数(
$ {N}_{e} $ )、RMSE随着粒子总数的增加而表现得更好,但处理多路径延迟量(M)时所需的时间也随之增加. 因此,适当增加粒子总数能够提高重采样粒子滤波的性能,有效去除伪距多路径序列中的噪声,但同时需要关注其对计算时间和资源的影响. 这一结论对于在实际应用中优化滤波算法参数具有重要的指导意义.4) 定位实验表明,本文提出的基于粒子滤波重采样的伪距多路径处理方法,可实现PPP定位精度E、N与U方向提升达10.68%、11.95%与26.13%的效果,动态定位精度E、N与U方向分别提升5.9%、1.1%与53.4%,对完善数据预处理模型具有重要作用.
因此,本文提出的修正方法不仅能够在一定程度上提升GNSS位置服务的性能,还有望为解决多路径误差问题提供新的思路和途径. 这对于提升GNSS系统的整体性能和可靠性具有重要意义.
致谢:感谢MGEX数据产品支持.
-
表 1 各系统伪距多路径误差统计
m 卫星系统 RMSE AVE 原始 改正 原始 改正 BDS 0.361 5 0.147 6 0.010 4 0.005 1 GPS 0.609 6 0.256 0 −0.262 0 0.195 0 GLONASS 0.430 9 0.165 1 0.004 8 0.001 9 Galileo 0.526 2 0.159 3 −0.010 8 0.010 1 表 2 各实测点伪距多路径误差统计
m 实测点 RMSE AVE 原始数据 改正数据 原始数据 改正数据 Hzjy 0.334 1 0.231 0 0.208 3 −0.151 7 Deqi 0.276 5 0.181 6 −0.108 3 0.088 2 Lina 0.291 3 0.183 5 −0.085 0 0.072 4 Fuya 0.273 2 0.191 5 0.115 8 −0.106 9 Xish 0.239 7 0.168 0 0.016 5 0.007 0 Keqi 0.500 2 0.178 9 0.002 1 −0.001 4 表 3 基于不同多路径处理模型的定位精度统计
方向 GPS GPS+BDS WN 粒子滤波 提升率/% WN 粒子滤波 提升率/% E 2.88 2.88 0.00 2.34 2.09 10.68 N 2.26 1.99 11.95 1.05 1.05 0.00 U 26.18 19.34 26.13 20.18 16.33 19.08 表 4 基于不同多路径处理模型的动态定位残差统计
多路径策略 E N U WN/m 0.051 0.020 0.641 粒子滤波/m 0.048 0.018 0.299 提升率/% 5.9 1.1 53.4 -
[1] ZHAN W, HE X F, JIA D Z, et al. Mitigation of multipath effects in GPS and BDS positioning using window matching method based sidereal filtering[J]. GPS solutions, 2022(57): 427-446. DOI: 10.1007/s40328-022-00384-6
[2] LIU Z F, TIAN Y M, XIONG W H, et al. A local filtering approach to mitigating the GNSS multipath effects in relative precise positioning considering the multipath spatial correlation[J]. Advances in space research, 2024, 74(6): 2709-2727. DOI: 10.1016/j.asr.2024.03.017
[3] ZHANG Z T. Code and phase multipath mitigation by using the observation-domain parameterization and its application in five-frequency GNSS ambiguity resolution[J]. GPS solution, 2021(25): 144. DOI: 10.1007/s10291-021-01179-y
[4] PAN Y X, MOLLER G, SOJA B. Machine learning-based multipath modeling in spatial domain applied to GNSS short baseline processing[J]. GPS solution, 2024, 28(1): 9. DOI: 10.1007/s10291-023-01553-y
[5] ZHANG Z T, YUAN H J, LI B F, et al. Feasibility of easy-to-implement methods to analyze systematic errors of multipath, differential code bias, and inter-system bias for low-cost receivers[J]. GPS solution, 2021, 25(3): 1-14. DOI: 10.1007/s10291-021-01149-4
[6] TIAN Y M, LIU Z F, LIN M, et al. Modelling and mitigation of GNSS multipath effects by least-squares collocation considering spatial autocorrelation[J]. Journal of geodesy, 2023(97): 37. DOI: 10.1007/s00190-023-01726-0
[7] WANNINGER L, BEER S. BeiDou satellite-induced code pseudorange variations: diagnosis and therapy[J]. GPS solution, 2015(19): 639-648. DOI: 10.1007/s10291-014-0423-3
[8] 王颖喆, 陶贤露, 朱锋, 等. 利用智能手机实现GNSS原始观测值的高精度差分定位[J]. 武汉大学学报(信息科学版), 2021, 46(12): 1941-1950. [9] 胡超, 王潜心, 郭忠臣, 等. 一种基于GNSS全系统全频点观测的多路径修正及定位模型[J]. 武汉大学学报(信息科学版), 2023: 1-13. [10] 龚学文, 王甫红. 星载GPS伪距多路径误差与观测噪声对自主定轨的影响分析[J]. 武汉大学学报(信息科学版), 2018, 43(7): 1048-1055. [11] 吴涛, 胡艳霞, 田甜, 等. GNSS干扰定位技术分析[J]. 全球定位系统, 2023, 48(5): 103-111. DOI: 10.12265/j.gnss.2023100 [12] ZENG J W, ZHANG Z T, HE X F, et al. Real-time GNSS multiple cycle slip detection and repair based on a controllable geometry-based method in relative positioning[J]. Measurement, 2023(216): 112940. DOI: 10.1016/j.measurement.2023.112940
[13] 张芮, 熊永良, 雷飞. 一种融合轨检信息的GNSS动态单历元多路径误差提取方法[J]. 武汉大学学报(信息科学版), 2021, 46(6): 905-912. [14] LU R, CHEN W, LI Z, et al. An improved joint modeling method for multipath mitigation of GPS, BDS-3, and Galileo overlapping frequency signals in typical environments[J]. Journal of geodesy, 2023, 97(10): 95. DOI: 10.57757/IUGG23-3840
[15] QIU W Q, ZENG Q H, XU R, et al. A multipath mitigation algorithm for GNSS signals based on the steepest descent approach[J]. Satellite navigation, 2022, 3(14): 80-90. DOI: 10.1007/s00190-023-01788-0
[16] LI X R, WANG L, QU X Y, et al. A GPS multipath mitigation method in coordinate-domain considering the effects of gross errors and missing data[J]. Measurement, 2024(225): 114035. DOI: 10.1016/j.measurement.2023.114035
[17] LU R, CHEN W, DONG D N, et al. Multipath mitigation in GNSS precise point positioning based on trend-surface analysis and multipath hemispherical map[J]. GPS solution, 2021, 25(3):14. DOI: 10.1007/s10291-021-01156-5
[18] LI X Z, XIONG Y L, XU S G, et al. A multipath error reduction method for BDS using Tikhonov regularization with parameter optimization[J]. Remote sensing, 2023, 15(13): 3400. DOI: 10.3390/rs15133400
[19] 周仁宇, 胡志刚, 蔡洪亮, 等. 使用抛物面定向天线分析北斗三号星上伪距和载波测距偏差[J]. 武汉大学学报(信息科学版), 2021, 46(9): 1298-1308. [20] SUZUKI T, MATSUO K, AMANO Y. Rotating GNSS antennas: simultaneous LOS and NLOS multipath mitigation[J]. GPS solution, 2020, 24(3): 86. DOI: 10.1007/s10291-020-01006-w
[21] PETER S, GARCIA-ASENJO L, SERGIO B. Optimal combination and reference functions of signal-to-noise measurements for GNSS multipath detection[J]. Measurement science and technology, 2019, 30(4):1-13. DOI: 10.1088/1361-6501/ab05ae
[22] 何聪聪, 王中元, 张坦. 北斗观测数据质量及动态PPP性能分析[J]. 全球定位系统, 2024, 49(1): 120-126. DOI: 10.12265/j.gnss.2023183 [23] ZHANG Z T, LI B F, GAO Y, et al. Real-time carrier phase multipath detection based on dual-frequency C/N0 data[J]. GPS solution, 2019, 23(7): 1-13. DOI: 10.1007/s10291-018-0799
[24] ZOU X, LI Z Y, WANG Y W, et al. Multipath error fusion modeling methods for multi-GNSS[J]. Remote sensing, 2021, 13(15): 2925. DOI: 10.3390/rs13152925
[25] FUHRMANN T, LUO X G, KNOPFLER A, et al. Generating statistically robust multipath stacking maps using congruent cells[J]. GPS solution, 2015, 19(1): 83-92. DOI: 10.1007/s10291-014-0367-7
[26] ZHANG X, ZHANG B C, YUAN Y B, et al. Extending multipath hemispherical model to account for time-varying receiver code biases[J]. Advances in space research, 2020(65): 650-662. DOI: 10.1016/j.asr.2019.11.003
[27] 张智超, 贾小林, 焦文海, 等. 抗差Vondrak滤波方法在时间频率传递中的应用[J]. 全球定位系统, 2024, 29(2):23-29. [28] ZHENG D W, ZHONG P, DING X L, et al. Filtering GPS time-series using a Vondrak filter and cross-validation[J]. Journal of geodesy, 2005(79): 363-369. DOI: 10.1007/s00190-005-0474-x
[29] 王笑蕾, 何秀凤, 陈殊, 等. 地基GNSS-IR风速反演原理及方法初探[J]. 测绘学报, 2021, 50(10): 1298-1307. DOI: 10.11947/j.AGCS.2021.20200586 [30] LIN X, LI W, LI S D, et al. Combined adaptive robust Kalman filter algorithm[J]. Measurement science and technology, 2021, 32(7). DOI: 10.1088/1361-6501/abf57c
[31] 王宪, 袁洪. 基于粒子滤波的GPS多径估计[J]. 控制与决策, 2010, 25(8): 1139-1143, 1148. [32] 宫轶松, 归庆明, 李保利, 等. 基于均值漂移的粒子滤波算法设计及其在导航数据处理中的应用[J]. 测绘学报, 2011(40): 120-125, 132. [33] 种衍文, 王泽文, 陈蓉, 等. 一种多特征自适应融合的粒子滤波红外目标跟踪方法[J]. 武汉大学学报(信息科学版), 2016, 14(5): 598-640. [34] 褚鼎立, 陈红, 王旭光. 基于自适应权重和模拟退火的鲸鱼优化算法[J]. 电子学报, 2019, 47(5): 992-999.