CN106324629A - 一种bds_gps_glonass融合精密单点定位方法 - Google Patents

Abstract

本发明提供了一种BDS_GPS_GLONASS融合精密单点定位方法，首先读取GPS、GLONASS、BDS三系统的原始数据，采用M‑W组合探测法进行周跳检测及粗差剔除，然后对三系统时空基准进行统一；接着构建三系统观测方程，采用双频伪距和载波相位观测值的消电离层组合方法进行误差改正，建立随机模型；最终采用扩展卡尔曼滤波进行模糊度和坐标求解。本发明可以大大改善卫星的几何图形强度，显著提高定位系统的自主完备性，有效缩短定位收敛时间以及提高在单系统卫星数较少或者卫星星座分布较差时的定位精度。

Description

一种BDS_GPS_GLONASS融合精密单点定位方法

技术领域

本发明涉及一种定位方法，属于卫星导航高精度非差定位领域。

背景技术

在传统GPS静态定位中，绝对定位的精度在20m左右，满足不了精密导航和测量用户的需求，一般都采用相对定位方式，通过组成双差观测值消除接收机钟差、卫星钟差等公共误差及削弱对流层延迟、电离层延迟等相关性强的误差影响，来达到提高精度的目的。这种观测方式无需考虑复杂的误差模型，具有解算模型简单、定位精度高等优势；但也存在一些不足，如观测时至少有一台接收机置于已知站上，影响了作业效率，提高了作业成本；另外，随着距离的增加，对流层延迟、电离层延迟等误差的相关性减弱，必须相应地延长观测时间，才能达到预期精度。精密单点定位技术，与单点定位相比，其静态定位精度可达到厘米级，定位精度高；与RTK技术相比，无需与其它组配或架设基准站，只需要单机作业，从而节约成本，且提高工作效率，同时不受距离约束，数据解算可快速完成，数据处理更加简单；其厘米级的定位结果又可以满足一般导航和测量用户的需求，因此，在测量的各领域有广泛的应用前景。当前的精密单点定位技术的实现主要是基于GPS观测数据。但GPS作为一种基于卫星的定位技术，系统的可用性、定位结果的可靠性和精度很大程度上取决于观测到的卫星的数量。随着北斗卫星导航系统(BDS)于2012年底正式提供亚太地区区域性服务，BDS加入GNSS行列，研究BDS与其他系统之间的组合定位尤为重要，特别是在单系统受所处环境影响导致可视卫星颗数较少时，组合系统能够有效的提高单系统的定位精度。利用多系统融合，可以增加观测卫星数，有助于提高定位精度，增加恶劣环境下的可用性。GPS+GLONASS+BDS在可用性、连续性和完好性方面的保障远比单一系统好，研究利用GPS+GLONASS+BDS观测数据进行精密单点定位，将使可视卫星数比单一GPS系统大大增加。首先，将有助于提高精密定位的精度；其次将有效的解决在城市、山区等遮挡严重的地方，观测GPS卫星数少于4颗，导航定位无法实现的问题，弥补了单一GPS系统卫星数目不足的问题；有效的降低用户对某一系统的依赖性，解决单一系统导航定位安全性不足的问题。

发明内容

为了克服现有技术的不足，本发明提供一种基于GPS+GLONASS+BDS融合的精密单点定位方法，三系统组合精密定位的时间系统和空间系统，建立了精密单点定位非差组合模型，相比单一系统而言，GPS+GLONASS+BDS无论是从连续性、可用性、可靠性、精度以及效率等各方面都更具优势。GPS、G LONASS及BDS作为三个完全独立的全球卫星导航定位系统，三者的融合可以大大改善卫星的几何图形强度，能显著提高定位系统的自主完备性，有效缩短定位收敛时间以及提高在单系统卫星数较少或者卫星星座分布较差时的定位精度。

本发明解决其技术问题所采用的技术方案包括以下步骤：

(1)读取GPS、GLONASS、BDS三系统的原始数据，包括RINEX格式原始观测数据、广播星历，精密星历和精密钟差；

(2)对GPS、GLONASS、BDS三系统的原始数据采用M-W组合探测法进行周跳检测及粗差剔除；用检测过周跳的载波相位来平滑剔除粗差的伪距，对GPS、GLONASS、BDS三系统时空基准进行统一处理；

(3)通过拉格朗日插值方法内插精密星历、钟差，构建GPS、GLONASS、BDS三系统观测方程，采用双频伪距和载波相位观测值的消电离层组合方法进行误差改正；

(4)对于同一系统的不同卫星，基于高度角的先验定权方法，建立随机模型；对于不同系统下的卫星，采用的是Helmert验后方差估计的方法进行定权，建立随机模型；

(5)采用扩展卡尔曼滤波进行模糊度和坐标求解。

本发明的有益效果是：融合定位系统收敛时间相比单系统收敛时间有所减少，这对于改善单系统精密单点定位需要长时间收敛时间具有重要意义。融合定位系统能够有效的提高单系统的可视卫星颗数，改善组合系统的空间卫星分布。特别是在环境恶劣情况下，单系统卫星颗数较少时，融合定位系统能够发挥明显优势。利用全天原始观测数据进行实施，全天时间中GPS平均可视卫星颗数5.7颗；BDS平均可视卫星数6.3颗；GLONASS平均可视卫星数5.0颗；三系统融合定位的平均可视卫星颗数最高，达到17颗。本发明能显著提高定位系统的自主完备性和精度。

附图说明

图1是本发明的流程示意图；

图2是实施例的定位结果X分量误差示意图；

图3是实施例的定位结果Y分量误差示意图；

图4是实施例的定位结果Z分量误差示意图。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

本发明包括以下步骤：

(1)读取原始数据

说明：原始数据包括：GPS、GLONASS、BDS三系统RINEX格式原始观测数据、广播星历，精密星历、精密钟差。

(2)数据预处理

说明：对GPS、GLONASS、BDS三系统原始数据采用M-W组合探测法进行周跳检测及粗差剔除。用检测过周跳的载波相位来平滑剔除粗差的伪距，对GPS、GLONASS、BDS三系统时空基准进行统一处理。

(3)建立精密单点定位模型

说明：通过拉格朗日插值方法内插精密星历、钟差，构建GPS、GLONASS、BDS三系统观测方程，采用双频伪距和载波相位观测值的消电离层组合方法进行误差改正。

(4)建立随机模型

说明：对于同一系统的不同卫星，基于高度角的先验定权方法，对于不同系统下的卫星，采用的是Helmert验后方差估计的方法来进行定权。

(5)解算结果

说明：采用扩展卡尔曼滤波(EKF)进行模糊度和坐标求解。

本发明的实施例如图1所示，其步骤如下：

步骤1：读取RINEX格式原始观测数据和广播星历文件，精密星历、精密钟差。

步骤2：对步骤1读取的数据进行数据预处理

(1)剔除粗差、周跳检测、平滑伪距

对所述原始观测数据依次进行剔除粗差、周跳检测、和平滑伪距处理，生成预处理后的原始观测数据。采用M-W组合探测法进行周跳检测及粗差剔除。用检测过周跳的载波相位来平滑剔除粗差的伪距，提高其精度，减少初始化的时间。

1)时空基准统一

√时间基准统一

GPS的时间系统(t_GPS)属于原子时系统，它的秒长和原子时秒长相同，但它和国际原子时存在着不同的原点，并且在任一时刻它们二者之间都存在一个常量偏差(19s)。北斗系统的时间基准为北斗时(t_BDS)。t_SDB采用国际单位制(SI)秒为基本单位连续累计，不闰秒，起始历元为2006年1月1日协调世界时(t_UTC)00时00分00秒，采用周和周内秒计数。t_BDS与t_UTC之间的闰秒信息在导航电文中传播。GLONASS时间系统(t_GLONASS)属于UTC时间系统，它是基于GLONASS同步中心(CS)时间产生的，由于GLONASS控制部分本身存在的特性，使得t_GLONASS与俄罗斯维持的协调世界时t_UTC(SU)存在的整数差为3h，此外它们之间还存在有1ms以内的系统误差τ_r。t_GPS与t_GLONASS之间的转换关系为：

t_GPS＝t_GLONASS+τ_c+τ_u+τ_g+τ_r

式中，τ_c＝t_UTC(SU)-t_GLONASS；τ_u＝t_UTC-t_UTC(SU)；τ_g＝t_GPS-t_UTCτ_r为GPS观测值与GLONASS观测值间的接收机钟差。

t_GPS与t_BDS之间的转换关系为：

t_GPS＝t_BDS+τ_c+τ_g+τ_r，τ_c＝t_UTC-t_BDS，τ_g＝t_GPS-t_UTC，τ_r为GPS观测值与BDS观测值间的接收机钟差。

√坐标基准统一

GPS坐标系统采用的是WGS-84坐标系。该坐标系是地心空间直角坐标系，它的原点为地球质心，Z轴指向国际时间服务机构(BIH)1984.0定义的协议地球极(CTP)方向，X轴指向(BIH)1984.0的零子午面和CTP赤道的交点，Y轴与Z轴、X轴垂直构成右手坐标系。BDS坐标系统采用的是我国2000大地坐标系统(CGCS2000)，它的定义和国际地球参考系统(ITRS)相一致。GLONASS的苏联1985标系(SGS85)，1993年以后改为使用PZ90.02坐标系，该坐标系是由俄罗斯进行地面网与空间网联合平差后所建立的。

GPS与BDS两个坐标系的主要不同在于参考椭球的扁率f有微小的差异，而这种差异在当前的测量精度水平中可以忽略，对于在坐标系定义上的比较，可以认为在同一历元下CGCS2000和WGS-84在坐标系的实现精度范围内，两者的坐标是一致的。当前，大多数学者认为使用俄罗斯MCC计算得到的七参数作为二个坐标系转换参数的精度最高，该七参数是通过全球激光跟踪测轨数据计算得到，PZ90.02与WGS-84之间的坐标转换公式为

${\left[\begin{array}{c}{X}_w\\ Y_w\\ Z_w\end{array}\right]}_{WGS-84}=\left[\begin{array}{c}-0.47\\ -0.51\\ -1.56\end{array}\right]+\left(1+22\×{10}^{-9}\right)\left[\begin{array}{ccc}1& -1.728\×{10}^{-6}& -0.017\×{10}^{-6}\\ 1.728\×{10}^{-6}& 1& 0.076\×{10}^{-6}\\ 0.017\×{10}^{-6}& -0.076\×{10}^{-6}& 1\end{array}\right]{\left[\begin{array}{c}{X}_p\\ Y_P\\ Z_P\end{array}\right]}_{PZ90.02}$

假设空间一点P在一个直角坐标系PZ90.02中的位置坐标为(x，y，z)，那么该点在另一个不同却几乎平行的直角坐标系WGS-84中的位置坐标(x′，y′，z′)为

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]=\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\\ \mathrm{\Δ}z\end{array}\right]+(1+{\mathrm{\δ}}_{PZWGS})\left[\begin{array}{ccc}1& {\mathrm{\δ}}_{\mathrm{\ω}}& -{\mathrm{\δ}}_{\mathrm{\ψ}}\\ -{\mathrm{\δ}}_{\mathrm{\ω}}& 1& {\mathrm{\δ}}_{\mathrm{\ϵ}}\\ {\mathrm{\δ}}_{\mathrm{\ψ}}& -{\mathrm{\δ}}_{\mathrm{\ϵ}}& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]$

坐标平移量(Δx，Δy，Δz)为坐标系PZ90.02的原点在坐标系WGS-84中的坐标，δ_ε、δ_ψ、δ_ω为坐标系PZ90.02分别绕X、Y、Z坐标轴旋转而能得到WGS-84相应

坐标轴的所需旋转角度，δ_PZWGS体现为坐标尺度因子。

步骤3：利用拉格朗日插值法内插精密星历和钟差，采用非差模型进行高精度单点定位，建立精密单点定位的误差方程如下：

$\begin{array}{c}{P}^g\left(L_i\right)={\mathrm{\ρ}}^g+c({\mathrm{dt}}^g-{\mathrm{dT}}^g)+{d^g}_{ord}+{d^g}_{trop}-{d^g}_{ion/Li}\\ +{d^g}_{mpath/P\left(Li\right)}+{\mathrm{\ϵ}}^g\left(P\left(L_i\right)\right)\end{array}---\left(1\right)$

P^l(L_i)＝ρ^l+c(dt^l-dT^l)+d^l _orb+d^l _trop-d^l _ion/Li+d^l _mpath/P(Li)+ε^l(P(L_i)) (3)

P^c(L_i)＝ρ^c+c(dt^c-dT^c)+d^c _orb+d^c _trop-d^c _ion/Li+d^c _mpath/P(Li)+ε^c(P(L_i)) (5)

Φ(L_i)是Li上的载波观测值(m)；ρ是几何距离，(X_s,Y_s,Z_s)是卫星发射时刻t的坐标，(x,y,z)是接收机坐标；c是光速(m/s)；dt是卫星钟差(s)；dT是接收机钟差(s)；d_orb是卫星轨道误差(m)；d_trop是对流层延迟(m)；d_ion/Li为Li上的电离层延迟(m)；λ_i是Li的波长(m)；N_i是Li的整周模糊度(周)；是接收机振荡器的初始相位；是卫星振荡器的初始相位；d_mpath/P(Li)是Li上的伪距测量值的多路径效应(m)；d_mpath/φ(Li)是Li上的载波相位测量值的多路径效应(m)；ε(·)是测量噪声。其中，g表示GPS卫星，L表示GLONASS卫星，c表示BDS卫星。

将对流层延迟d_trop表示成对流层延迟zpd与其映射函数M的积。采用双频无电离层非差观测值组合来消除电离层的一阶影响，定位的数学模型将变为：

$\begin{array}{c}{P}_{IF}^g=\frac{f_1^2{P}^g\left(L_1\right)-f_2^2{P}^g\left(L_2\right)}{f_1^2-f_2^2}={\mathrm{\ρ}}^g-{\mathrm{cdT}}^g+{d^g}_{ord}\\ +M\·zpd+d_{mpath/P^g\left(L1+L2\right)}+\mathrm{\ϵ}\left(P^g\left(L1+L2\right)\right)\end{array}---\left(7\right)$

$\begin{array}{c}{\mathrm{\Φ}}_{IF}^g=\frac{f_1^2{\mathrm{\φ}}^g\left(L_1\right)-f_2^2{\mathrm{\φ}}^g\left(L_2\right)}{f_1^2-f_2^2}={\mathrm{\ρ}}^g-{\mathrm{cdT}}^g+M\·zpd+{\mathrm{\λ}}^g{N}^g\\ +d_{mpath/{\mathrm{\φ}}^g\left(L1+L2\right)}+\mathrm{\ϵ}\left({\mathrm{\φ}}^g\left(L1+L2\right)\right)\end{array}---\left(8\right)$

$\begin{array}{c}{P}_{IF}^1=\frac{f_1^2{P}^l\left(L_1\right)-f_2^2{P}^l\left(L_2\right)}{f_1^2-f_2^2}={\mathrm{\ρ}}^l-{\mathrm{cdT}}^l+{d^l}_{orb}\\ +M\·zpd+d_{mpath/P^l\left(L1+L2\right)}+\mathrm{\ϵ}\left(P^l\left(L1+L2\right)\right)\end{array}---\left(9\right)$

$\begin{array}{c}{\mathrm{\Φ}}_{IF}^l=\frac{f_1^2{\mathrm{\φ}}^l\left(L_1\right)-f_2^2{\mathrm{\φ}}^l\left(L_2\right)}{f_1^2-f_2^2}={\mathrm{\ρ}}^l-{\mathrm{cdT}}^l+M\·zpd+{\mathrm{\λ}}^l{N}^l\\ +d_{mpath/{\mathrm{\φ}}^l\left(L1+L2\right)}+\mathrm{\ϵ}\left({\mathrm{\φ}}^l\left(L1+L2\right)\right)\end{array}---\left(10\right)$

$\begin{array}{c}{P}_{IF}^c=\frac{f_1^2{P}^c\left(L_1\right)-f_2^2{P}^c\left(L_2\right)}{f_1^2-f_2^2}={\mathrm{\ρ}}^c-{\mathrm{cdT}}^c+{d^c}_{orb}\\ +M\·zpd+d_{mpath/P^c\left(L1+L2\right)}+\mathrm{\ϵ}\left(P^c\left(L1+L2\right)\right)\end{array}---\left(11\right)$

$\begin{array}{c}{\mathrm{\Φ}}_{IF}^c=\frac{f_1^2{\mathrm{\φ}}^c\left(L_1\right)-f_2^2{\mathrm{\φ}}^c\left(L_2\right)}{f_1^2-f_2^2}={\mathrm{\ρ}}^c-{\mathrm{cdT}}^c+M\·zpd+{\mathrm{\λ}}^c{N}^c\\ +d_{mpath/{\mathrm{\φ}}^c\left(L1+L2\right)}+\mathrm{\ϵ}\left({\mathrm{\φ}}^c\left(L1+L2\right)\right)\end{array}---\left(12\right)$

步骤4：进行定位解算时，必须对伪距测量值进行电离层、相对论效应等误差改正，从而消除绝大部分伪距测量值和载波相位观测值的误差，但有部分误差有所残留，鉴于不同系统的不同伪距误差，就有必要对每颗卫星的权值进行设定。对于同一系统的不同卫星，基于高度角的先验定权方法，对于不同系统下的卫星，采用的是Helmert验后方差估计的方法来进行定权。观测值的精度可以通过高度角来间接反映。当出现卫星高度角偏低的情况，电离层、对流层、多路径等误差就对观测值出现比较大的影响，这就出现精度降低的情况。所以必须适当地提高卫星高度角。以卫星高度角定权的一种随机模型形式为：

σ² _C＝Sigma²c/sin²(E)，

σ² _Φ＝Sigma²Φ/sin²(E)

式中：下标C表示为伪距观测方程；下标Φ表示为相位观测方程；E为卫星高度角；Sigma表示观测值天顶方向测量精度。计算时观测值先验精度设置为：GPS：SigmaC：3m，SigmaΦ：0.02m；GLONASS：Sigma C：3.5m；SigmaΦ：0.025m；BDS：Sigma C：3.2m，SigmaΦ：0.02m。

步骤5：三系统融合定位中，观测方程为非线性，采用扩展卡尔曼滤波(EKF)进行模糊度和坐标积累，快速收敛。构造方程如下：

系统方程：

x_i＝(x_i,y_i,z_i,zpd_i,N_i) (13)

测量方程：

y_i＝(P_IF,i,Φ_IF,i) (14)

系统预测方程：

${\hat{x}}_i(+)={\hat{x}}_i(-)+K_i(y_i-h({\hat{x}}_i(-)\left)\right)---\left(15\right)$

系统预测方程协方差：

$P_i(+)=(I-K_{i}H({\hat{x}}_i(-)\left)\right)P_i(-)---\left(16\right)$

增益方程：

$K_i=P_i(-)H\left(\hat{x}\right(-\left)\right){\left(H\right({\hat{x}}_i(-)\left)P_i\right(-)H{\left({\hat{x}}_i(-)\right)}^T+R_i)}^{-1}---\left(17\right)$

式中：与P_i分别为在历元t_k时的待估计状态量及其相关系数矩阵；(—)与(+)分别为EKF观测更新前与观测更新后；h(x)、H(x)与R(x)分别为观测模型量、偏导数和观测误差的协方差矩阵。

当观测方程经过线性化后，状态量与其相关系数矩阵的时间更新为

$\left\{\begin{array}{c}{\hat{x}}_{i+1}(-)=F_i^{i+1}{\hat{x}}_i(+)\\ P_{i+1}(-)=F_i^{i+1}{P}_i(+)F_i^{i+1}+Q_i^{i+1}\end{array}\right.$

式中，与分别为历元t_k到t_k+1时系统噪声的转移矩阵及其相关系数矩阵。

定位结果如图2、3、4所示：融合定位系统能够减少收敛时间；融合系统的定位精度比单系统的定位精度高。

Claims (1)

1.一种BDS_GPS_GLONASS融合精密单点定位方法，其特征在于包括下述步骤：

(1)读取GPS、GLONASS、BDS三系统的原始数据，包括RINEX格式原始观测数据、广播星历，精密星历和精密钟差；

(2)对GPS、GLONASS、BDS三系统的原始数据采用M-W组合探测法进行周跳检测及粗差剔除；用检测过周跳的载波相位来平滑剔除粗差的伪距，对GPS、GLONASS、BDS三系统时空基准进行统一处理；

(3)通过拉格朗日插值方法内插精密星历、钟差，构建GPS、GLONASS、BDS三系统观测方程，采用双频伪距和载波相位观测值的消电离层组合方法进行误差改正；

(4)对于同一系统的不同卫星，基于高度角的先验定权方法，建立随机模型；对于不同系统下的卫星，采用的是Helmert验后方差估计的方法进行定权，建立随机模型；

(5)采用扩展卡尔曼滤波进行模糊度和坐标求解。