CN104597465A - 一种提高gps与glonass组合精密单点定位收敛速度的方法 - Google Patents

Abstract

本发明公开了一种提高GPS与GLONASS组合精密单点定位收敛速度的方法，通过在函数模型中引入多个独立的“时频偏差”参数，弥补了现有技术在建模中忽略GLONASS码频间偏差的缺陷，能在不损失定位精度的前提下，明显提高GPS与GLONASS组合PPP的收敛速度。采用本发明所提供的方法，在进行组合PPP时对码伪距观测值可采用GLONASS与GPS等权的随机模型，避免了现有技术依赖经验确定权比的不足。

Description

一种提高GPS与GLONASS组合精密单点定位收敛速度的方法

技术领域

本发明涉及卫星大地测量及导航定位技术领域，尤其涉及一种提高GPS与GLONASS组合精密单点定位收敛速度的方法。

背景技术

随着GLONASS现代化进程的不断推进，截至2011年12月8日，GLONASS在轨工作卫星达24颗，已恢复系统完全运行能力。与GPS不同，当前的GLONASS系统使用频分多址(Frequency Division Multiple Access,FDMA)技术区分来自不同卫星的信号。卫星信号在接收机内部频率通道会产生相应的硬件延迟，一般将其分为平均延迟偏差项和因通道频率而异的延迟偏差项，后者又称为频间偏差(Inter-Frequency Bias,IFB)。由于GLONASS采用FDMA技术，不同频率的卫星信号会进入接收机内部对应的频率通道，由此会产生不同的载波相位和码频间偏差。精密单点定位(Precise Point Positioning,PPP)技术在不固定整周模糊度的情况下，载波相位频间偏差将由模糊度参数吸收，因而无需在观测模型中考虑。

由于GPS与GLONASS存在系统时间基准的差异，现有GPS/GLONASS组合PPP技术主要采用“单参数”法构建函数模型，即在观测方程中仅估计GPS接收机钟差，同时引入一个“系统时差”(inter-system bias,ISB)参数来表征GPS和GLONASS系统时间偏差。现有技术的主要缺陷在于：引入的“系统时差”参数仅吸收接收机码平均硬件延迟，而GLONASS码频间偏差因各卫星频率间的差异而无法被接收机钟差参数吸收，最终将体现在码伪距观测值的残差中。众所周知，在不固定整周模糊度的情况下，初始阶段PPP浮点解模糊度的收敛速度主要取决于消电离层码伪距组合观测值的精度。组合PPP模型中忽略码频间偏差的影响将造成GLONASS码伪距观测值精度要低于其实际观测精度，尤其当码频间偏差在数值上远大于码伪距观测值自身的噪声水平时，GPS/GLONASS组合系统PPP的收敛速度甚至不如GPS单系统PPP。现有技术一般采用在随机模型中对GLONASS伪距观测值进行降权的方法以削弱上述影响，其局限性在于：要确定GPS与GLONASS码伪距观测值的合理权比，需要对测站接收机的GLONASS码频间偏差特性有一定了解，而目前绝大多数GNSS接收机生产商并不提供该参数的先验值，在实际操作中多依赖于数据处理人员的经验来确定权比，因而具有一定的盲目性，理论上不严密。

发明内容

本发明所要解决的技术问题是针对背景技术中所涉及的缺陷，提供一种提高GPS与GLONASS组合精密单点定位收敛速度的方法。

本发明为解决上述技术问题采用以下技术方案：

一种提高GPS与GLONASS组合精密单点定位收敛速度的方法，包括以下步骤：

步骤1)，获取GPS、GLONASS观测数据；

步骤2)，对步骤1)中获取的GPS、GLONASS观测数据进行周跳探测，并标记发生周跳的历元；

步骤3)，对步骤1)获取的GPS、GLONASS观测数据分别组成消电离层组合观测值；

步骤4)，分别建立GPS、GLONASS卫星观测方程，其中将GLONASS接收机码频间偏差参数与系统时差参数进行合并，为每颗观测的GLONASS卫星均设置一个独立时频偏差参数；

步骤5)，获取精密卫星轨道及卫星钟差数据，对观测方程中的各项误差源进行建模改正；

步骤6)，将所有GPS、GLONASS卫星观测方程线性化，并表示成矩阵形式；

步骤7)，采用扩展卡尔曼滤波方法进行参数估计；其中，对标记为发生周跳的历元，将模糊度参数重置后再进行扩展卡尔曼滤波估计；

步骤8)，输出组合PPP结果。

作为本发明一种提高GPS与GLONASS组合精密单点定位收敛速度的方法进一步的优化方案，步骤3)中对于GPS卫星i组成消电离层组合观测值的计算公式为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},i}^G=\frac{{f_1}^2{P_1}^G-{f_2}^2{P_2}^G}{{f_1}^2-{f_2}^{22}}\\ {\mathrm{\Φ}}_{\mathrm{IF},i}^G=\frac{{f_1}^2{{\mathrm{\Φ}}_1}^G-{f_2}^2{{\mathrm{\Φ}}_2}^G}{{f_1}^2-{f_2}^2}\end{array}\right.$

式中，G表示GPS系统，分别为码伪距和载波相位的消电离层组合观测值；f₁、f₂分别表示L₁和L₂载波的频率；P₁ ^G、P₂ ^G与Φ₁ ^G、Φ₂ ^G分别为相应频率的码伪距和载波相位原始观测值；

对于GLONASS卫星j组成消电离层组合观测值的计算公式为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},j}^R=\frac{f_{K,1}^2{P_1}^R-f_{K,2}^2{P_2}^R}{f_{K,1}^2-f_{K,2}^2}\\ {\mathrm{\Φ}}_{\mathrm{IF},j}^R=\frac{f_{K,1}^2{{\mathrm{\Φ}}_1}^R-f_{K,2}^2{{\mathrm{\Φ}}_2}^R}{f_{K,1}^2-f_{K,2}^2}\end{array}\right.$

式中，R表示GLONASS系统，分别为码伪距和载波相位的消电离层组合观测值；K表示卫星j所对应的GLONASS频率通道号；f_K,1、f_K,2分别表示卫星j所对应的GLONASS频率通道K的L₁和L₂载波的频率；P₁ ^R、P₂ ^R与Φ₁ ^R、Φ₂ ^R分别为测站r相应频率的码伪距和载波相位原始观测值。

作为本发明一种提高GPS与GLONASS组合精密单点定位收敛速度的方法进一步的优化方案，步骤4)中对于GPS卫星i的观测方程表示为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},i}^G={\mathrm{\ρ}}^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\mathrm{\δ}}_{{\mathrm{trop}}_i}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\ϵ}}_P^G\\ {\mathrm{\Φ}}_{\mathrm{IF},i}^G={\mathrm{\ρ}}^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\mathrm{\λ}}^G{N}_i^G+{\mathrm{\δ}}_{{\mathrm{trop}}_i}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\δ}}_{\mathrm{phw}}+{\mathrm{\ϵ}}_{\mathrm{\Φ}}^G\end{array}\right.$

式中，ρ^G为GPS卫地距；为待估的GPS接收机钟差参数；为精密卫星钟差；为GPS卫星i对应的消电离层相位模糊度；δ_tropi为对流层斜延迟；δ_tide为潮汐效应改正；δ_rel为相对论效应改正；δ_phw为天线相位缠绕改正；为观测噪声；c为真空中光速，λ^G为GPS消电离层相位观测值所对应的波长；

对于GLONASS卫星j，其观测方程表示为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},j}^R={\mathrm{\ρ}}^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\mathrm{\β}}_j+{\mathrm{\δ}}_{{\mathrm{trop}}_j}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\ϵ}}_P^R\\ {\mathrm{\Φ}}_{\mathrm{IF},j}^R={\mathrm{\ρ}}^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\mathrm{\β}}_j+{\mathrm{\λ}}_j^R{N}_j^R+{\mathrm{\δ}}_{{\mathrm{trop}}_j}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\ϵ}}_{\mathrm{\Φ}}^R\end{array}\right.$

式中，为GLONASS精密卫星钟差；为GLONASS卫星j对应的消电离层相位模糊度；为卫星j消电离层相位观测值所对应的波长；与GPS观测方程符号相同的参数其含义是相同的，此处用上标R表示GLONASS系统，β_j为待估的GLONASS卫星j的时频偏差ISFB参数，具体表示为：

${\mathrm{\β}}_j={\mathrm{dt}}_{\mathrm{SYS}}+M_r^R-D_r^G+{\stackrel{\‾}{B}}_{j,r}^R,(j=\mathrm{1,2},...,n)$

式中，dt_SYS为系统时间偏差；为GLONASS接收机码平均延迟偏差；为GPS接收机码延迟偏差；为GLONASS接收机码频间偏差，具体表示为：

${\stackrel{\‾}{B}}_{j,r}^R=\frac{f_{K,1}^2{\hat{B}}_{j,r}^{P1}-f_{K,2}^2{\hat{B}}_{j,r}^{P2}}{f_{K,1}^2-f_{K,2}^2}$

式中，为原始观测值对应的GLONASS接收机码频间偏差。

作为本发明一种提高GPS与GLONASS组合精密单点定位收敛速度的方法进一步的优化方案，步骤6)中所有卫星的观测方程线性化后的矩阵形式为：

$\left\{\begin{array}{c}y\left(t\right)=\mathrm{AX}\left(t\right)+{\mathrm{\ϵ}}_y\\ {\mathrm{\ϵ}}_y~N(0,{\mathrm{\Ω}}_y)\end{array}\right.$

式中，y(t)为观测向量，其维数为2(m+n)；m、n分别为在历元t、测站r同步观测到的GPS卫星、GLONASS卫星的数量；X(t)为待估参数向量，其维数为(m+2n+5)，包括测站位置[X Y Z]^T、GPS接收机钟差天顶对流层湿延迟δ_zwd、GPS相位模糊度 ${\left[\begin{array}{ccc}{N}_1^G& ...& N_m^G\end{array}\right]}^T,$ GLONASS相位模糊度 ${\left[\begin{array}{ccc}{N}_1^R& ...& N_n^R\end{array}\right]}^T$ 及时频偏差[β₁ … β_n]^T；A为设计阵；ε_y为测量噪声向量，且服从均值为0、协方差矩阵为Ω_y的正态分布；Ω_y为对角矩阵，其对角线元素表示为卫星高度角的三角函数，由下式计算：

${\mathrm{\σ}}^2=\frac{{\mathrm{\σ}}_0^2}{{\sin }^{2}E}$

式中，σ为观测值先验中误差；E为测站处的卫星高度角；σ₀为天顶方向观测值的先验中误差。

作为本发明一种提高GPS与GLONASS组合精密单点定位收敛速度的方法进一步的优化方案，步骤7)中所述扩展卡尔曼滤波的状态方程表示如下：

X(t)＝Ψ_t,t-1X(t-1)+w_t-1,w_t-1～(0,Ω_w)

式中，X(t)、X(t-1)为相邻历元的系统状态向量；Ψ_t,t-1为状态转移矩阵；w_t-1为状态噪声向量，且服从均值为0、协方差矩阵为Ω_w的正态分布；

状态转移矩阵Ψ_t,t-1表示为：

${\mathrm{\Ψ}}_{t,t-1}=\left[\begin{array}{ccccc}\underset{3\×3}{P}& & & & \\ & 0& & & \\ & & 1& & \\ & & & \underset{(m+n)\×(m+n)}{I}& \\ & & & & \underset{n\×n}{I}\end{array}\right]$

式中，P为位置参数转移矩阵，I为单位阵；

协方差矩阵Ω_w表示为：

${\mathrm{\Ω}}_w=\left[\begin{array}{ccccc}\underset{3\×3}{q_{\mathrm{pos}}}\mathrm{\Δt}& & & & \\ & q_{{\mathrm{dt}}_r^G}\mathrm{\Δt}& & & \\ & & q_{{\mathrm{\δ}}_{\mathrm{zwd}}}\mathrm{\Δt}& & \\ & & & \underset{(m+n)\×(m+n)}{O}& \\ & & & & \underset{n\×n}{q_{\mathrm{\β}}}\mathrm{\Δt}\end{array}\right]$

式中，Δt为相邻历元的时间间隔，q_pos、q_β分别表示待求测站位置、GPS接收机钟差、对流层天顶湿延迟及时频偏差参数的谱密度矩阵，O表示零矩阵。

作为本发明一种提高GPS与GLONASS组合精密单点定位收敛速度的方法进一步的优化方案，所述步骤4)中GLONASS观测方程中接收机钟差参数与GPS观测方程中的设置相同，且引入的时频偏差参数由GPS与GLONASS的系统时间偏差、GLONASS接收机码平均延迟、需扣除的GPS接收机码延迟及GLONASS接收机码频间偏差四部分组成。

作为本发明一种提高GPS与GLONASS组合精密单点定位收敛速度的方法进一步的优化方案，所述GPS可见卫星数m≥5。

作为本发明一种提高GPS与GLONASS组合精密单点定位收敛速度的方法进一步的优化方案，将时频偏差参数变化视为随机游走过程处理。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

与现有技术相比，本发明通过在函数模型中引入多个独立的“时频偏差”(ISFB)参数，弥补了现有技术在建模中忽略GLONASS码频间偏差的不足，能在不损失定位精度的前提下，明显提高GPS/GLONASS组合PPP的收敛速度；采用本发明所提供的方法，实现了对GLONASS码频间偏差的函数模型补偿，因此在进行组合PPP时对码伪距观测值采用GLONASS与GPS等权的随机模型亦能获得良好效果，避免了现有技术依赖经验确定权比的盲目性。

附图说明

图1是本发明方法的流程图；

图2是三种实验方案3D位置收敛时间对比图；

图3是实验方案2和方案3的GLONASS码残差分布对比图；

图4所有样本测站三种实验方案的收敛时间对比图；

图5所有样本测站三种实验方案的定位结果对比图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明，但是本发明的保护范围不局限于所述实施例。

如图1所示，本发明提供的一种提高GPS/GLONASS组合精密单点定位收敛速度的方法，具体步骤如下：

步骤1：获取GPS、GLONASS观测数据；为检验本发明所提供方法的有效性，采用来自欧洲区域IGS参考站网和加拿大NRCan网供12个测站数据进行实验；所有数据均由GPS/GLONASS双系统接收机采集，共涉及6个GNSS设备商的9种不同接收机类型，具体信息如表1所示。数据采集时间为2013年年积日(Day of Year,DOY)91～180，采样率为30s，观测时长为24h。

表1：实验所采用配备相同品牌接收机的测站信息

步骤2：对步骤1中获取的GPS/GLONASS观测数据进行周跳探测，并标记发生周跳的历元；

步骤3：对步骤1获取的GPS/GLONASS观测数据分别组成消电离层组合观测值，通过以下公式计算得到：

对于GPS卫星i：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},i}^G=\frac{{f_1}^2{P_1}^G-{f_2}^2{P_2}^G}{{f_1}^2-{f_2}^{22}}\\ {\mathrm{\Φ}}_{\mathrm{IF},i}^G=\frac{{f_1}^2{{\mathrm{\Φ}}_1}^G-{f_2}^2{{\mathrm{\Φ}}_2}^G}{{f_1}^2-{f_2}^2}\end{array}\right.---\left(1\right),$

式中，分别为码伪距和载波相位的消电离层组合观测值；f₁、f₂分别表示L₁和L₂载波的频率；P₁ ^G、P₂ ^G与Φ₁ ^G、Φ₂ ^G分别为相应频率的码伪距和载波相位原始观测值；

对于GLONASS卫星j：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},j}^R=\frac{f_{K,1}^2{P_1}^R-f_{K,2}^2{P_2}^R}{f_{K,1}^2-f_{K,2}^2}\\ {\mathrm{\Φ}}_{\mathrm{IF},j}^R=\frac{f_{K,1}^2{{\mathrm{\Φ}}_1}^R-f_{K,2}^2{{\mathrm{\Φ}}_2}^R}{f_{K,1}^2-f_{K,2}^2}\end{array}\right.---\left(2\right),$

式中，与GPS计算公式符号相同的参数其含义是相同的，此处用上标R表示GLONASS系统，下标K表示卫星j所对应的GLONASS频率通道号；

步骤4：分别建立GPS、GLONASS卫星观测方程，其中将GLONASS接收机码“频间偏差”与“系统时差”参数进行合并，为每颗观测的GLONASS卫星均设置一个独立“时频偏差”(Inter-System and Inter-Frequency Bias,ISFB)参数，具体表示如下：

对于GPS卫星i，其观测方程表示为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},i}^G={\mathrm{\ρ}}^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\mathrm{\δ}}_{{\mathrm{trop}}_i}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\ϵ}}_P^G\\ {\mathrm{\Φ}}_{\mathrm{IF},i}^G={\mathrm{\ρ}}^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\mathrm{\λ}}^G{N}_i^G+{\mathrm{\δ}}_{{\mathrm{trop}}_i}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\δ}}_{\mathrm{phw}}+{\mathrm{\ϵ}}_{\mathrm{\Φ}}^G\end{array}\right.---\left(3\right),$

式中，ρ^G为GPS卫地距；为待估的GPS接收机钟差参数；为精密卫星钟差；为GPS卫星i对应的消电离层相位模糊度；为对流层斜延迟；δ_tide为潮汐效应改正；δ_rel为相对论效应改正；δ_phw为天线相位缠绕改正；为观测噪声；

对于GLONASS卫星j，其观测方程表示为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},j}^R={\mathrm{\ρ}}^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\mathrm{\β}}_j+{\mathrm{\δ}}_{{\mathrm{trop}}_j}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{+\mathrm{\ϵ}}_P^R\\ {\mathrm{\Φ}}_{\mathrm{IF},j}^R={\mathrm{\ρ}}^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\mathrm{\β}}_j+{\mathrm{\λ}}_j^R{N}_j^R+{\mathrm{\δ}}_{{\mathrm{trop}}_j}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\ϵ}}_{\mathrm{\Φ}}^R\end{array}\right.---\left(4\right),$

式中，与GPS观测方程符号相同的参数其含义是相同的，此处用上标R表示GLONASS系统，β_j为待估的GLONASS卫星j的时频偏差ISFB参数，具体表示为：

${\mathrm{\β}}_j={\mathrm{dt}}_{\mathrm{SYS}}+M_r^R-D_r^G+{\stackrel{\‾}{B}}_{j,r}^R,(j=\mathrm{1,2},...,n)---\left(5\right),$

式中，dt_SYS为系统时间偏差；为GLONASS接收机码平均延迟偏差；为GPS接收机码延迟偏差；为GLONASS接收机码频间偏差，具体表示为：

${\stackrel{\‾}{B}}_{j,r}^R=\frac{f_{K,1}^2{\hat{B}}_{j,r}^{P1}-f_{K,2}^2{\hat{B}}_{j,r}^{P2}}{f_{K,1}^2-f_{K,2}^2}---\left(6\right),$

式中，为原始观测值对应的GLONASS接收机码频间偏差，f_K,1、f_K,2为相应的信号频率，对应的GLONASS频率通道号为K；

注意本步骤中GLONASS观测方程中接收机钟差参数与GPS观测方程中的相同，而GPS与GLONASS系统时间差及GLONASS码频间偏差均由时频偏差(ISFB)参数吸收，并作为待估参数进行估计。

步骤5：获取精密卫星轨道及卫星钟差数据，对观测方程中的各项误差源进行建模改正；GPS和GLONASS精密卫星轨道及卫星钟差均统一采用欧洲空间局(European space agency,ESA)提供的事后产品，其中卫星轨道与钟差产品采样间隔分别为15min、30s；对流层干延迟采用Saastamoinen模型及GMF映射函数进行改正；地球固体潮、海洋负荷潮及极潮参照IERS Conventions(2010)提供的模型计算，其中海潮模型采用FES2004；相对论效应改正参照GPS系统接口规范IS-GPS-200G提供的模型计算；地球自转参数(ERP)采用IGS提供的事后产品；不同类型码观测量的卫星硬件延迟偏差(DCB)采用欧洲定轨中心(CODE)提供的事后产品；卫星与接收机相位天线偏差及变化改正采用igs08.atx文件中提供数值，并且进行天线相位缠绕的改正。

步骤6：将所有GPS/GLONASS卫星观测方程线性化，并表示成矩阵形式；假设在历元t，测站r同步观测到m颗GPS卫星、n颗GLONASS卫星，将所有卫星的观测方程线性化，写成矩阵形式表示为：

$\left\{\begin{array}{c}y\left(t\right)=\mathrm{AX}\left(t\right)+{\mathrm{\ϵ}}_y\\ {\mathrm{\ϵ}}_y~N(0,{\mathrm{\Ω}}_y)\end{array}\right.---\left(7\right),$

式中，y(t)为观测向量，其维数为2(m+n)；X(t)为待估参数向量，其维数为(m+2n+5)，包括测站位置[X Y Z]^T、GPS接收机钟差天顶对流层湿延迟δ_zwd、GPS相位模糊度 ${\left[\begin{array}{ccc}{\tilde{N}}_1^G& ...& {\tilde{N}}_m^G\end{array}\right]}^T,$ GLONASS相位模糊度 ${\left[\begin{array}{ccc}{\tilde{N}}_1^R& ...& {\tilde{N}}_n^R\end{array}\right]}^T$ 及时频偏差[β₁ … β_n]^T；A为设计阵；要获得逐个历元的定位解，需GPS可见卫星数m≥5；ε_y为测量噪声向量，且服从均值为0、协方差矩阵为Ω_y的正态分布；Ω_y为对角矩阵，其对角线元素表示为卫星高度角的三角函数，由下式计算：

${\mathrm{\σ}}^2=\frac{{\mathrm{\σ}}_0^2}{{\sin }^{2}E}---\left(8\right),$

式中，σ为观测值先验中误差；E为测站处的卫星高度角；σ₀为天顶方向观测值的先验中误差。

步骤7：采用扩展卡尔曼滤波(Extended Kalman Filter,EKF)模型进行参数估计，其中对标记为发生周跳历元的模糊度参数进行重置。EKF状态方程表示如下：

X(t)＝Ψ_t,t-1X(t-1)+w_t-1,w_t-1～(0,Ω_w)        (9)，

式中，X(t)、X(t-1)为相邻历元的系统状态向量；Ψ_t,t-1为状态转移矩阵；w_t-1为状态噪声向量，且服从均值为0、协方差矩阵为Ω_w的正态分布。

状态转移矩阵Ψ_t,t-1表示为：

${\mathrm{\Ψ}}_{t,t-1}=\left[\begin{array}{ccccc}\underset{3\×3}{P}& & & & \\ & 0& & & \\ & & 1& & \\ & & & \underset{(m+n)\×(m+n)}{I}& \\ & & & & \underset{n\×n}{I}\end{array}\right]---\left(10\right),$

式中，P为位置参数转移矩阵，I为单位阵。

协方差矩阵Ω_w表示为：

${\mathrm{\Ω}}_w=\left[\begin{array}{ccccc}\underset{3\×3}{q_{\mathrm{pos}}}\mathrm{\Δt}& & & & \\ & q_{{\mathrm{dt}}_r^G}\mathrm{\Δt}& & & \\ & & q_{{\mathrm{\δ}}_{\mathrm{zwd}}}\mathrm{\Δt}& & \\ & & & \underset{(m+n)\×(m+n)}{O}& \\ & & & & \underset{n\×n}{q_{\mathrm{\β}}}\mathrm{\Δt}\end{array}\right]---\left(11\right),$

式中，Δt为相邻历元的时间间隔，q表示各状态参数的谱密度矩阵。在随机模型中，将GPS接收机钟差视为白噪声处理，其先验中误差设置为10²m；天顶对流层湿延迟视为随机游走过程处理，其先验中误差设置为由于相邻历元间隔Δt通常较小，在此较短时间内可以认为ISFB参数变化是稳定的，因而同样将其视为随机游走过程处理，其先验中误差设置为在处理静态数据时，接收机位置参数视为常量；GPS与GLONASS模糊度参数均视为不随时间变化的量。

为分析本发明提供PPP方法的效果，设计3种方案对实验数据进行计算。

方案1：GPS单系统PPP；

方案2：采用传统“单参数”法的GPS/GLONASS组合PPP，其中GPS与GLONASS载波相位观测值权比按等权处理，而对GLONASS码伪距观测值进行降权处理，二者的先验中误差之比为1:1.5；

方案3：采用本发明“多参数”法的GPS/GLONASS组合PPP，其中GPS与GLONASS载波相位观测值先验中误差设置与方案2相同，对码伪距观测值采用GLONASS与GPS等权的随机模型。

需要说明的是，对于PPP所涉及到各项误差源如卫星位置误差、卫星钟差、潮汐效应、相对论效应、相位缠绕、天线相位偏差及变化等，3种方案的改正策略均设置为相同。

图2为2013年DOY 91测站LAMA、ZIM2及KUUJ的3D位置偏差的收敛情况。从图2可以看出，对于测站LAMA、ZIM2，方案2的收敛速度明显优于方案1，这说明GLONASS观测值的加入对PPP有明显的收敛加速效果；而对于测站KUUJ，方案2的收敛时间明显大于方案1，这意味着GLONASS观测值的加入反而不如仅使用GPS单系统观测值收敛速度快；方案3通过设置多个“时频偏差”参数吸收各GLONASS卫星对应的码频间偏差，函数模型中GLONASS码伪距观测值实际精度可以认为与其自身的观测精度已非常接近，此时设置GPS与GLONASS观测值等权亦能取得较好效果。显然，图2中所有3个测站方案3的收敛时间都是最短的，其中对测站KUUJ收敛的加速效果尤为明显。

图3进一步给上述3个测站所有GLONASS卫星的码伪距残差统计分布。从图3可以看出，方案3所得到的GLONASS码残差分布较方案2有显著改善，各测站码残差均值在0附近，且标准差约为方案2的1/2左右，符合较为理想的正态分布；这说明本发明“多参数”法对于GLONASS码频间偏差的函数模型补偿是有效的。

为进一步验证本发明“多参数”法对提高GPS/GLONASS组合PPP收敛速度的有效性，对实验数据涉及的12个测站均采用上述3种方案对2013年DOY 91～180的观测数据进行处理；每3h时段数据为一组统计各方案的收敛时间，收敛水平分别设置为10cm和5cm，结果如图4所示。从图4可以看出，在3D偏差收敛至10cm水平，除配备TPS品牌接收机的测站方案1和方案2的收敛时间相当外，其他5种品牌方案2相对于方案1的收敛时间可缩短30％左右；采用方案3后，所有测站收敛时间进一步明显缩短，与方案2相比，配备TPS品牌接收机测站ALGO、KUUJ及配备JPS品牌接收机测站ONSA的收敛速度均提高近50％，其他测站提高约20％；在3D偏差收敛至5cm水平，采用方案3后，测站ALGO、KUUJ、ONSA的收敛速度相对于方案2还可进一步提高约20％。

本发明方法引入多个附加参数将使函数模型的自由度减小，为分析是否会对组合PPP的定位精度带来不利影响，对实验数据按3种方案分别处理得到3h和24h定位结果，定位结果3D偏差RMS值如图5所示。从图5可以看出，对于3h时段，方案2、方案3的定位精度相当，二者要明显优于方案1，仅有数个测站方案3的RMS值要略大于方案2；对于24h时段，3种方案在不同测站的表现各有优劣，整体上RMS值互差为毫米级，因而可以认为3种方案单天解精度是相当的。由此可见，方案3虽引入多个附加参数将使函数模型的自由度减小，但实际并未明显降低组合PPP的定位精度，与方案2的RMS值差异在毫米级水平。

本发明的具体实施方式中未涉及的说明属于本领域公知的技术，可参考公知技术加以实施。

Claims (8)

1.一种提高GPS与GLONASS组合精密单点定位收敛速度的方法，其特征在于，包括以下步骤：

步骤1)，获取GPS、GLONASS观测数据；

步骤2)，对步骤1)中获取的GPS、GLONASS观测数据进行周跳探测，并标记发生周跳的历元；

步骤3)，对步骤1)获取的GPS、GLONASS观测数据分别组成消电离层组合观测值；

步骤4)，分别建立GPS、GLONASS卫星观测方程，其中将GLONASS接收机码频间偏差参数与系统时差参数进行合并，为每颗观测的GLONASS卫星均设置一个独立时频偏差参数；

步骤5)，获取精密卫星轨道及卫星钟差数据，对观测方程中的各项误差源进行建模改正；

步骤6)，将所有GPS、GLONASS卫星观测方程线性化，并表示成矩阵形式；

步骤7)，采用扩展卡尔曼滤波方法进行参数估计；其中，对标记为发生周跳的历元，将模糊度参数重置后再进行扩展卡尔曼滤波估计；

步骤8)，输出组合PPP结果。

2.根据权利要求1所述的提高GPS与GLONASS组合精密单点定位收敛速度的方法，其特征在于，步骤3)中对于GPS卫星i组成消电离层组合观测值的计算公式为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},i}^G=\frac{f_1^2{P}_1^G-f_2^2{P_2}^G}{f_1^2-f_2^2}\\ {\mathrm{\Φ}}_{\mathrm{IF},i}^G=\frac{f_1^2{\mathrm{\Φ}}_1^G-f_2^2{\mathrm{\Φ}}_2^G}{f_1^2-f_2^2}\end{array}\right.$

式中，G表示GPS系统，分别为码伪距和载波相位的消电离层组合观测值；f₁、f₂分别表示L₁和L₂载波的频率；与分别为相应频率的码伪距和载波相位原始观测值；

对于GLONASS卫星j组成消电离层组合观测值的计算公式为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},j}^R=\frac{f_{K,1}^2{P}_1^R-f_{K,2}^2{P_2}^R}{f_{K,1}^2-f_{K,2}^2}\\ {\mathrm{\Φ}}_{\mathrm{IF},j}^R=\frac{f_{K,1}^2{\mathrm{\Φ}}_1^R-f_{K,2}^2{\mathrm{\Φ}}_2^R}{f_{K,1}^2-f_{K,2}^2}\end{array}\right.$

式中，R表示GLONASS系统，分别为码伪距和载波相位的消电离层组合观测值；K表示卫星j所对应的GLONASS频率通道号；f_K,1、f_K,2分别表示卫星j所对应的GLONASS频率通道K的L₁和L₂载波的频率；与分别为测站r相应频率的码伪距和载波相位原始观测值。

3.根据权利要求2所述的提高GPS与GLONASS组合精密单点定位收敛速度的方法，其特征在于，步骤4)中对于GPS卫星i的观测方程表示为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},i}^G={\mathrm{\ρ}}^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\mathrm{\δ}}_{{\mathrm{trop}}_i}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\ϵ}}_P^G\\ {\mathrm{\Φ}}_{\mathrm{IF},i}^G={\mathrm{\ρ}}^G+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_i^{G/s})+{\mathrm{\λ}}^G{N}_i^G+{\mathrm{\δ}}_{{\mathrm{trop}}_i}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\δ}}_{\mathrm{phw}}+{\mathrm{\ϵ}}_{\mathrm{\Φ}}^G\end{array}\right.$

式中，ρ^G为GPS卫地距；为待估的GPS接收机钟差参数；为GPS精密卫星钟差；为GPS卫星i对应的消电离层相位模糊度；为对流层斜延迟；δ_tide为潮汐效应改正；δ_rel为相对论效应改正；δ_phw为天线相位缠绕改正；为观测噪声；c为真空中光速，λ^G为GPS消电离层相位观测值所对应的波长；

对于GLONASS卫星j，其观测方程表示为：

$\left\{\begin{array}{c}{P}_{\mathrm{IF},j}^R={\mathrm{\ρ}}^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\mathrm{\β}}_j+{\mathrm{\δ}}_{{\mathrm{trop}}_j}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\ϵ}}_P^R\\ {\mathrm{\Φ}}_{\mathrm{IF},j}^R={\mathrm{\ρ}}^R+c({\mathrm{dt}}_r^G-{\mathrm{dt}}_j^{R/s})+{\mathrm{\β}}_j+{\mathrm{\λ}}_j^R{N}_j^G+{\mathrm{\δ}}_{{\mathrm{trop}}_j}+{\mathrm{\δ}}_{\mathrm{tide}}+{\mathrm{\δ}}_{\mathrm{rel}}+{\mathrm{\δ}}_{\mathrm{phw}}+{\mathrm{\ϵ}}_{\mathrm{\Φ}}^R\end{array}\right.$

式中，为GLONASS精密卫星钟差；为GLONASS卫星j对应的消电离层相位模糊度；为卫星j消电离层相位观测值所对应的波长；与GPS观测方程符号相同的参数其含义是相同的，此处用上标R表示GLONASS系统，β_j为待估的GLONASS卫星j的时频偏差ISFB参数，具体表示为：

${\mathrm{\β}}_j={\mathrm{dt}}_{\mathrm{SYS}}+M_r^R-D_r^G+{\stackrel{\‾}{B}}_{j,r}^R,(j=\mathrm{1,2},...,n)$

式中，dt_SYS为系统时间偏差；为GLONASS接收机码平均延迟偏差；为GPS接收机码延迟偏差；为GLONASS接收机码频间偏差，具体表示为：

${\stackrel{\‾}{B}}_{j,r}^R=\frac{f_{K,1}^2{\hat{B}}_{j,r}^{P1}-f_{K,2}^2{\hat{B}}_{j,r}^{P2}}{f_{K,1}^2-f_{K,2}^2}$

式中，为原始观测值对应的GLONASS接收机码频间偏差。

4.根据权利要求3所述的提高GPS与GLONASS组合精密单点定位收敛速度的方法，其特征在于，步骤6)中所有卫星的观测方程线性化后的矩阵形式为：

$\left\{\begin{array}{c}y\left(t\right)=\mathrm{AX}\left(t\right)+{\mathrm{\ϵ}}_y\\ {\mathrm{\ϵ}}_y~N(0,{\mathrm{\Ω}}_y)\end{array}\right.$

式中，y(t)为观测向量，其维数为2(m+n)；m、n分别为在历元t、测站r同步观测到的GPS卫星、GLONASS卫星的数量；X(t)为待估参数向量，其维数为(m+2n+5)，包括测站位置[X Y Z]^T、GPS接收机钟差天顶对流层湿延迟δ_zwd、GPS相位模糊度GLONASS相位模糊度及时频偏差[β₁ … β_n]^T；A为设计阵；ε_y为测量噪声向量，且服从均值为0、协方差矩阵为Ω_y的正态分布；Ω_y为对角矩阵，其对角线元素表示为卫星高度角的三角函数，由下式计算：

${\mathrm{\σ}}^2=\frac{{\mathrm{\σ}}_0^2}{{\sin }^{2}E}$

式中，σ为观测值先验中误差；E为测站处的卫星高度角；σ₀为天顶方向观测值的先验中误差。

5.根据权利要求4所述的提高GPS与GLONASS组合精密单点定位收敛速度的方法，其特征在于，步骤7)中所述扩展卡尔曼滤波的状态方程表示如下：

X(t)＝Ψ_t,t-1X(t-1)+w_t-1,w_t-1～(0,Ω_w)

式中，X(t)、X(t-1)为相邻历元的系统状态向量；Ψ_t,t-1为状态转移矩阵；w_t-1为状态噪声向量，且服从均值为0、协方差矩阵为Ω_w的正态分布；

状态转移矩阵Ψ_t,t-1表示为：

${\mathrm{\Ψ}}_{t,t-1}=\left[\begin{array}{ccccc}\underset{3\×3}{P}& & & & \\ & 0& & & \\ & & 1& & \\ & & & \underset{(m+n)\×(m+n)}{I}& \\ & & & & \underset{n\×n}{I}\end{array}\right]$

式中，P为位置参数转移矩阵，I为单位阵；

协方差矩阵Ω_w表示为：

${\mathrm{\Ω}}_w=\left[\begin{array}{ccccc}\underset{3\×3}{q_{\mathrm{pos}}\mathrm{\Δt}}& & & & \\ & q_{{\mathrm{dt}}_r^G}\mathrm{\Δt}& & & \\ & & q_{{\mathrm{\δ}}_{\mathrm{zwd}}}\mathrm{\Δt}& & \\ & & & \underset{(m+n)\×(m+n)}{O}& \\ & & & & \underset{n\×n}{q_{\mathrm{\β}}}\mathrm{\Δt}\end{array}\right]$

式中，Δt为相邻历元的时间间隔，q_pos、q_β分别表示待求测站位置、GPS接收机钟差、对流层天顶湿延迟及时频偏差参数的谱密度矩阵，O表示零矩阵。

6.根据权利要求1所述的提高GPS与GLONASS组合精密单点定位收敛速度的方法，其特征在于，所述步骤4)中GLONASS观测方程中接收机钟差参数与GPS观测方程中的设置相同，且引入的时频偏差参数由GPS与GLONASS的系统时间偏差、GLONASS接收机码平均延迟、需扣除的GPS接收机码延迟及GLONASS接收机码频间偏差四部分组成。

7.根据权利要求4所述的提高GPS与GLONASS组合精密单点定位收敛速度的方法，其特征在于，所述GPS可见卫星数m≥5。

8.根据权利要求5所述的一种提高GPS与GLONASS组合精密单点定位收敛速度的方法，其特征在于，将时频偏差参数变化视为随机游走过程处理。