CN105182388A - 一种快速收敛的精密单点定位方法 - Google Patents

Abstract

本发明公开了一种快速收敛的精密单点定位方法，首先对全球导航卫星的观测数据进行预处理，再联合原始载波相位观测值的观测方程、码相位半合组合观测值的观测方程和几何无关组合观测值的观测方程构建改进的精密单点定位的观测方程，采用卡尔曼滤波对改进的精密单点定位的观测方程中的参数进行估计，得到接收机位置信息和对流层延迟、电离层延迟信息，实现GNSS快速单点精密定位。与现有技术相比：本发明一种快速收敛的精密单点定位方法能更好地降低观测噪声、轨道误差对滤波收敛速度与定位精度影响，并利用对流层延迟、电离层延迟历元间缓慢变化的特性作为约束条件，提高了精密单点定位解算的性能。

Description

一种快速收敛的精密单点定位方法

技术领域

本发明涉及全球导航卫星定位技术领域，尤其涉及一种快速收敛的精密单点定位方法。

背景技术

精密单点定位(PrecisePointPositioning，简称PPP)是基于状态空间域改正信息的高精度定位方法，基本实现方法为：用户采用一台GNSS(GlobalNavigationSatelliteSystem)接收机的双频码伪距和载波相位观测值，通过国际GNSS服务组织(InternationalGNSSService，简称IGS)提供的GNSS精密轨道和精密钟差产品内插出相应观测时刻的卫星轨道和钟差信息，同时应用误差模型修正天线相位中心、地球自转、固体潮汐等误差源的影响，进行单站绝对定位。目前，位置解收敛到静态厘米级、动态分米级的精度通常需要较长的时间，如静态模式下一般需要20～30分钟才能收敛到厘米级。收敛速度慢是限制PPP技术更加广泛应用的主要原因。

自从精密单点定位技术被提出以来，许多专家学者提出了多种定位方法，如Kouba在期刊“GPSSolutions”2001年第5卷第2期中提出了标准消电离层组合方法(Un-Differenceionosphere-freecombinedmethod，简称UD方法)，GaoY在期刊“Navigation”2002年第49卷第2期中提出了UofC方法(UniversityofCalgarymethod，简称UofC方法)。这两种方法均采用了消电离层组合，虽然降低了电离层延迟对定位的影响，但放大了观测噪声和多路径效应，不利于位置解的快速收敛；张宝成在“测绘学报”2010年第39卷第5期中提出了非差非组合方法(Un-Combinedmethod，简称UC方法)，该方法采用参数估计形式将电离层延迟模拟成随机游走，但采用了原始码伪距观测值，观测噪声仍然较大，且无法降低轨道误差对定位的影响。另外，GeM在期刊“JournalofGeodesy”2008年82卷第7期中、张小红在期刊“武汉大学学报·信息科学版”2010年第35卷第6期中、ShiJ在期刊“GPSSolutions”2014年第66卷第1期中均实现了基于相位偏差改正的精密单点定位固定解，虽在一定程度上提高了PPP收敛速度和定位精度，但参数的解算仍需基于上述三种基本方法来实现。

发明内容

本发明所要解决的技术问题是提供一种快速收敛的精密单点定位方法，有效降低观测噪声、轨道误差对滤波收敛速度与定位精度影响，并利用对流层延迟、电离层延迟历元间缓慢变化的特性作为约束条件，提高了精密单点定位解算的性能。

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

一种快速收敛的精密单点定位方法,包括如下步骤:

(S1)、对全球导航卫星的观测数据进行预处理：

采用M-W组合法和电离层残差法对原始载波相位观测值中存在的周跳进行探测，并修复；

对天线相位中心偏差、地球自转、相对论效应、天线相位缠绕误差源的影响进行修正；

对精密卫星轨道和卫星钟差数据进行拉格朗日多项式内插，获取观测时刻的卫星轨道和卫星钟差；

(S2)、顾及电离层延迟对不同频率观测值的影响关系，建立卫星的原始载波相位观测值和的观测方程，如公式(5)、公式(6)所示：

${\mathrm{\Φ}}_1^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_1^s}---\left(5\right)$

${\mathrm{\Φ}}_2^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-f_1^2/f_2^2\·I_1^s+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_2^s}---\left(6\right)$

其中，s表示卫星，和表示载波相位观测值，ρ表示卫星与接收机之间的几何距离，c表示真空中光的速度，表示接收机钟差，表示卫星钟差，分别表示对流层干、湿迟延映射函数；δ_zhd、δ_zwd分别表示对流层天顶方向干延迟和湿延迟，为载波相位观测值中站星视线方向的电离层延迟，分别为原始载波相位观测值和的模糊度，f₁、f₂为载波相位的频率，λ₁、λ₂表示载波相位的波长，表示载波相位观测值的观测噪声；

(S3)、根据原始码伪距观测值和原始载波相位观测值构建码相位半合组合观测值的观测方程，如公式(7)所示：

$L_{P_2{\mathrm{\Φ}}_2}^s=(P_2^s+{\mathrm{\Φ}}_2^s)/2=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s.{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s/2+{\mathrm{\ϵ}}_{L_{P_2{\mathrm{\Φ}}_2}^s}---\left(7\right)$

(S4)、根据原始载波相位观测值和构建几何无关组合观测值的观测方程，如公式(8)所示：

${\mathrm{\Φ}}_{GF}^s={\mathrm{\Φ}}_1^s-{\mathrm{\Φ}}_2^s=(f_1^2/f_2^2-1)\·I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s-{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_{GF}^s}---\left(8\right)$

(S5)、联合原始载波相位观测值和的观测方程、码相位半合组合观测值的观测方程和几何无关组合观测值的观测方程构建改进的精密单点定位的观测方程，如公式(9)所示，

$\left\{\begin{array}{c}{\mathrm{\Φ}}_1^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_2^s}\\ {\mathrm{\Φ}}_2^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-f_1^2/f_2^2\·I_1^s+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_2^s}\\ L_{P_2{\mathrm{\Φ}}_2}^s=\left(P_2^s+{\mathrm{\Φ}}_2^s\right)/2=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s/2+{\mathrm{\ϵ}}_{L_{P_2{\mathrm{\Φ}}_2}^s}\\ {\mathrm{\Φ}}_{GF}^s={\mathrm{\Φ}}_1^s-{\mathrm{\Φ}}_2^s=\left(f_1^s/f_2^2-1\right)\·I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s-{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_{GF}^s}\end{array}\right.---\left(9\right)$

(S6)、采用卡尔曼滤波对改进的精密单点定位的观测方程中的参数进行估计，得到接收机位置信息和对流层延迟、电离层延迟信息。

作为上述方案的进一步优化，对步骤(S5)的改进的精密单点定位的观测方程进行线性化处理，得到公式(10)：

y(k)＝A·X(k)+ε_y，ε_y～N(0,Ω_y)(10)

其中，y(k)为观测向量，其维数为4j，j为历元k同步观测到的卫星数量；X(k)为待估参数向量，维数为(5+3j)，A为设计矩阵，ε_y为测量噪声向量，Ω_y为测量噪声向量ε_y的协方差阵；y(k)、X(k)和A表达式分别为

$y\left(k\right)={\left[\begin{array}{cccccccccccc}{\mathrm{\Φ}}_1^1& \mathrm{...}& {\mathrm{\Φ}}_1^j& {\mathrm{\Φ}}_2^1& \mathrm{...}& {\mathrm{\Φ}}_2^j& L_{P_2{\mathrm{\Φ}}_2}^1& \mathrm{...}& L_{P_2{\mathrm{\Φ}}_2}^j& {\mathrm{\Φ}}_{GF}^1& \mathrm{...}& {\mathrm{\Φ}}_{GF}^j\end{array}\right]}^T$

$X\left(k\right)={\left[\begin{array}{cccccccccccccc}X& Y& Z& c\·\mathrm{\δ}{\stackrel{\‾}{t}}_r& {\mathrm{\δ}}_{zpw}& I_1^1& \mathrm{...}& I_1^j& {\stackrel{\‾}{N}}_1^1& \mathrm{...}& {\stackrel{\‾}{N}}_1^j& {\stackrel{\‾}{N}}_2^1& \mathrm{...}& {\stackrel{\‾}{N}}_2^j\end{array}\right]}^T$

$A=\[e\begin{array}{ccc}{}_4\⊗B& C\⊗I_j& D\⊗I_j\end{array}\]$

其中，X、Y、Z为接收机三维位置，表示克罗内克积，e₄为各元素均为1的4维列向量，I_j为j维单位矩阵，B、C和D子矩阵的含义分别为:

$B=\left[\begin{array}{ccccc}{l}^1& m^1& n^1& M_w^1& 1\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ l^j& m^j& n^j& M_w^j& 1\end{array}\right],C=\left[\begin{array}{c}-1\\ -f_1^2/f_2^2\\ 0\\ f_1^2/f_2^2-1\end{array}\right],D=\left[\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\\ 0& {\mathrm{\λ}}_2/2\\ {\mathrm{\λ}}_1& -{\mathrm{\λ}}_2\end{array}\right]$

其中，l、m、n分别为接收机到卫星的方向余弦。

作为上述方案的进一步优化，步骤(S5)得到的改进的精密单点定位的观测方程第k历元对应的协方差阵Ω_y(k)，如公式(11)所示：

${\mathrm{\Ω}}_y\left(k\right)=F\⊗I_j---\left(11\right)$

其中， $F=\left[\begin{array}{cccc}{\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2& 0& 0& {\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2\\ 0& {\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2& \frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}& -{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2\\ 0& \frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}& \frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2+{\mathrm{\σ}}_{P_2}^2}{4}& -\frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}\\ {\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2& -{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2& -\frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}& {\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2+{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2\end{array}\right],$ 表示克罗内克积，I_j为j维单位矩阵，分别为卫星原始载波相位观测值Φ₁、Φ₂，原始码伪距观测值P₂的测量噪声的方差。

与已有技术相比，本发明的一种快速收敛的精密单点定位方法的有益效果体现在：

1、本发明的一种快速收敛的精密单点定位方法有效地降低观测噪声、轨道误差对滤波收敛速度与定位精度影响，提高了精密单点定位解算的性能。

2、本发明的一种快速收敛的精密单点定位方法，不仅降低了观测噪声和轨道误差对定位的影响，还可以利用大气延迟平稳变化的特性作为约束提高滤波收敛速度，其构建的改进的精密单点定位的观测方程，具有更好的模型结构，采用卡尔曼滤波对改进的精密单点定位的观测方程进行参数估计，能克服现有解算方法的技术缺陷。

3、本发明的一种快速收敛的精密单点定位方法，能快速地实现精密单点定位厘米级位置解，可提高精密单点定位收敛速度以及定位精度。

附图说明

图1是本发明的一种快速收敛的精密单点定位方法的流程图。

图2为本发明的应用实施例采用的IGS观测站的分布图。

图3为使用本发明方法和UC方法对公共滤波收敛数据处理后的收敛时间统计结果图。

图4为使用本发明方法和UC方法对公共滤波收敛数据处理后的平面位置偏差统计结果图。

图5为使用本发明方法和UofC方法对公共滤波收敛数据处理后的收敛时间统计结果图。

图6为使用本发明方法和UofC方法对公共滤波收敛数据处理后的平面位置偏差统计结果图。

图7为使用本发明方法和UD方法对公共滤波收敛数据处理后的收敛时间统计结果图。

图8为使用本发明方法和UD方法对公共滤波收敛数据处理后的平面位置偏差统计结果图。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

一种快速收敛的精密单点定位方法,包括如下步骤:参见图1，图1是本发明的一种快速收敛的精密单点定位方法的流程图。

(S1)、对全球导航卫星的观测数据进行预处理：

采用M-W组合法和电离层残差法对原始载波相位观测值中存在的周跳进行探测，并修复；

对天线相位中心偏差、地球自转、相对论效应、天线相位缠绕误差源的影响进行修正；

对精密卫星轨道和卫星钟差数据进行拉格朗日多项式内插，获取观测时刻的卫星轨道和卫星钟差。

(S2)、建立卫星的原始载波相位观测值和的观测方程：

在GNSS应用中，在卫星s的原始码伪距和载波相位观测值的观测方程表达式如下：

$P_i^s=\mathrm{\ρ}+c\·{\mathrm{\δt}}_r-c\·{\mathrm{\δt}}^s+M_d^s\·{\mathrm{\δ}}_{zhd}+M_w^s\·{\mathrm{\δ}}_{zwd}+I_i^s+c(d_{r,i}-d_i^s)+{\mathrm{\ϵ}}_{P_i^s}---\left(1\right)$

其中，s表示卫星(s＝1,…,j)，i表示载波相位和码伪距的频率数，表示码伪距观测值，表示以米为单位的载波相位观测值，ρ表示卫星与接收机之间的几何距离，c表示真空中光的速度，δt_r表示接收机钟差，δt^s表示卫星钟差，分别表示对流层干、湿迟延映射函数，δ_zhd、δ_zwd分别表示对流层天顶方向干延迟和湿延迟，表示站星视线方向的电离层延迟，d_r,i表示接收机端的硬件码延迟，表示卫星端的硬件码延迟，表示载波相位的模糊度，λ_i表示载波相位的波长，表示接收机端的相位延迟，表示卫星端的相位延迟，分别表示码伪距和载波相位观测值的观测噪声；

将式(1)和式(2)简化为

$P_i^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s.{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}+I_i^s+{\mathrm{\ϵ}}_{P_i^s}---\left(3\right)$

${\mathrm{\Φ}}_i^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s.{\mathrm{\δ}}_{zpw}-I_i^s+{\mathrm{\λ}}_i\·{\stackrel{\‾}{N}}_i^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_i^s}---\left(4\right)$

其中，分别表示重新定义后的接收机钟差、卫星钟差和载波相位模糊度，表达式分别为

顾及电离层延迟对不同频率观测值的影响关系为则得到原始载波相位观测值和的观测方程，表达式为

${\mathrm{\Φ}}_1^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_1^s}---\left(5\right)$

${\mathrm{\Φ}}_2^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-f_1^2/f_2^2\·I_1^s+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_2^s}---\left(6\right)$

式(5)或式(6)中，和表示载波相位观测值，f₁、f₂为载波相位的频率，λ₁、λ₂表示载波相位的波长，表示载波相位观测值的观测噪声；I₁为载波相位观测值中站星视线方向的电离层延迟，分别为原始载波相位观测值和的模糊度；

(S3)、根据原始码伪距观测值和原始载波相位观测值构建码相位半合组合观测值的观测方程，如公式(7)所示：

$L_{P_2{\mathrm{\Φ}}_2}^s=(P_2^s+{\mathrm{\Φ}}_2^s)/2=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s.{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s/2+{\mathrm{\ϵ}}_{L_{P_2{\mathrm{\Φ}}_2}^s}---\left(7\right)$

(S4)、根据原始载波相位观测值和构建几何无关组合观测值的观测方程，如公式(8)所示：

${\mathrm{\Φ}}_{GF}^s={\mathrm{\Φ}}_1^s-{\mathrm{\Φ}}_2^s=(f_1^2/f_2^2-1)\·I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s-{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_{GF}^s}---\left(8\right)$

(S5)、联合原始载波相位观测值和的观测方程、码相位半合组合观测值的观测方程和几何无关组合观测值的观测方程构建改进的精密单点定位的观测方程，如公式(9)所示：

$\left\{\begin{array}{c}{\mathrm{\Φ}}_1^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_2^s}\\ {\mathrm{\Φ}}_2^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-f_1^2/f_2^2\·I_1^s+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_2^s}\\ L_{P_2{\mathrm{\Φ}}_2}^s=\left(P_2^s+{\mathrm{\Φ}}_2^s\right)/2=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s/2+{\mathrm{\ϵ}}_{L_{P_2{\mathrm{\Φ}}_2}^s}\\ {\mathrm{\Φ}}_{GF}^s={\mathrm{\Φ}}_1^s-{\mathrm{\Φ}}_2^s=\left(f_1^s/f_2^2-1\right)\·I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s-{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_{GF}^s}\end{array}\right.---\left(9\right)$

式(9)线性化后的矩阵形式为：

y(k)＝A·X(k)+ε_y，ε_y～N(0,Ω_y)(10)

式(10)中，y(k)为观测向量，其维数为4j，j为历元k同步观测到的卫星数量；X(k)为待估参数向量，维数为(5+3j)，A为设计矩阵，ε_y为测量噪声向量，Ω_y为测量噪声向量ε_y的协方差阵；y(k)、X(k)和A表达式分别为

$y\left(k\right)={\left[\begin{array}{cccccccccccc}{\mathrm{\Φ}}_1^1& \mathrm{...}& {\mathrm{\Φ}}_1^j& {\mathrm{\Φ}}_2^1& \mathrm{...}& {\mathrm{\Φ}}_2^j& L_{P_2{\mathrm{\Φ}}_2}^1& \mathrm{...}& L_{P_2{\mathrm{\Φ}}_2}^j& {\mathrm{\Φ}}_{GF}^1& \mathrm{...}& {\mathrm{\Φ}}_{GF}^j\end{array}\right]}^T,$

$X\left(k\right)={\left[\begin{array}{cccccccccccccc}X& Y& Z& c\·\mathrm{\δ}{\stackrel{\‾}{t}}_r& {\mathrm{\δ}}_{zpw}& I_1^1& \mathrm{...}& I_1^j& {\stackrel{\‾}{N}}_1^1& \mathrm{...}& {\stackrel{\‾}{N}}_1^j& {\stackrel{\‾}{N}}_2^1& \mathrm{...}& {\stackrel{\‾}{N}}_2^j\end{array}\right]}^T,$

$A=\[e\begin{array}{ccc}{}_4\⊗B& C\⊗I_j& D\⊗I_j\end{array}\],$

其中，X、Y、Z为接收机三维位置，表示克罗内克积，e₄为各元素均为1的4维列向量，I_j为j维单位矩阵，B、C和D子矩阵的含义分别为：

$B=\left[\begin{array}{ccccc}{l}^1& m^1& n^1& M_w^1& 1\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ l^j& m^j& n^j& M_w^j& 1\end{array}\right],C=\left[\begin{array}{c}-1\\ -f_1^2/f_2^2\\ 0\\ f_1^2/f_2^2-1\end{array}\right],D=\left[\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\\ 0& {\mathrm{\λ}}_2/2\\ {\mathrm{\λ}}_1& -{\mathrm{\λ}}_2\end{array}\right],$

其中，l、m、n分别为接收机到卫星的方向余弦。

观测值协方差阵衡量了观测值的精度及相关性。由于卫星观测值受到噪声及多路径效应影响主要发生在低高度角卫星，为了不使高仰角卫星的观测值降低权重，精密单点定位中一般采用高度角定权。不考虑码和载波之间的相关性，即则步骤(S5)得到的改进的精密单点定位的观测方程第k历元对应的协方差阵Ω_y(k)，如公式(11)所示：

${\mathrm{\Ω}}_y\left(k\right)=F\⊗I_j---\left(11\right)$

式(11)中， $F=\left[\begin{array}{cccc}{\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2& 0& 0& {\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2\\ 0& {\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2& \frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}& -{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2\\ 0& \frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}& \frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2+{\mathrm{\σ}}_{P_2}^2}{4}& -\frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}\\ {\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2& -{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2& -\frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}& {\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2+{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2\end{array}\right],$ 表示克罗内克积，I_j为j维单位矩阵，分别为卫星原始载波相位观测值Φ₁、Φ₂，原始码伪距观测值P₂的测量噪声的方差。

(S6)、采用卡尔曼滤波对改进的精密单点定位的观测方程中的参数进行估计，得到接收机位置信息和对流层延迟、电离层延迟信息，实现GNSS快速单点精密定位。

在PPP数据处理中，通常采用卡尔曼(Kalman)滤波对参数进行估计，状态方程表达为

X_k+1＝Φ_k+1,kX_k+w_k,w_k～N(0,Q_k)(12)

式(12)中，X_k表示状态向量，Φ_k+1,k为状态转移矩阵，w_k为动态噪声，Q_k为动态噪声w_k的协方差阵。

动态噪声协方差阵Q_k反映了各待估参数的时变特性。在本发明方法中，主要参数包括接收机三维位置、接收机钟差、对流层天顶湿延迟、视线方向电离层延迟和载波相位模糊度,其中接收机钟差、对流层天顶湿延迟、视线方向电离层延迟模型化为随机游走，载波相位模糊度在参数估计时当作常数处理。于是,各待估参数动态噪声如下

$Q_k=\left[\begin{array}{ccccc}\underset{3\×3}{q_p}\·\mathrm{\Δ}t& 0& 0& 0& 0\\ 0& q_{t_r}\·\mathrm{\Δ}t& 0& 0& 0\\ 0& 0& q_{trop}\·\mathrm{\Δ}t& 0& 0\\ 0& 0& 0& \underset{n\×n}{q_{ion}}\·\mathrm{\Δ}t& 0\\ 0& 0& 0& 0& \underset{n\×n}{0}\end{array}\right]---\left(13\right)$

式(13)中,q_p为接收机三维位置谱密度,在静态和动态定位中分别为单位阵和零矩阵,为接收机钟差谱密度；q_trop为对流层天顶湿延迟的谱密度，q_ion为视线方向电离层延迟的谱密度，大气误差谱密度的选取与其变化率有关，Δt为相邻历元的时间间隔。

应用实施例

测试站点选取全球范围分布的144个IGS连续运行观测站，测站分布如图2，图2为本发明的应用实施例采用的IGS观测站的分布图。选取2013年6月17日00时至24时的观测数据，采样间隔为30s，间隔为15min的精密星历和间隔为5min的钟差产品，卫星截至高度角设置为5°。对144个IGS站24小时的观测数据以每0.5小时进行截取，得6912组数据。采用IGS发布的各个观测站坐标作为参考值，当水平位置偏差从某一历元开始均小于10cm时，认为滤波收敛，实现厘米级定位，而从第一个历元到该历元的时间间隔为收敛时间。为分析本发明提供PPP的定位方法的效果，处理步骤为：

步骤(1)：对应用实施例中的6912组数据分别使用本发明的方法、UC方法、UofC方法、UD方法进行计算，对平面位置收敛到厘米级精度的百分比和平均收敛时间进行分析；得到本发明的方法、UC方法、UofC方法和UD方法的收敛百分比和收敛时间统计信息，参见表1。

从表1中，可以得到：本发明方法具有最高的滤波收敛比例和最少的收敛时间，分别为70.2％和12.9min，较UC方法、UofC方法和UD方法收敛比例提高了19.2、2.8和49.3个百分点，平均收敛时间分别缩短了18.4％、7.2％和29.5％。可见，本发明方法能够较高比例的收敛到厘米级，且收敛时间较短。

表1收敛百分比和收敛时间统计信息

步骤(2)：选取本发明方法与UC方法公共滤波收敛部分，使用本发明方法与UC方法对收敛时间和平面位置偏差方面进行分析，得到收敛时间统计结果和平面位置偏差统计结果，参见图3和图4，图3为使用本发明方法和UC方法对公共滤波收敛数据处理后的收敛时间统计结果图；图4为使用本发明方法和UC方法对公共滤波收敛数据处理后的平面位置偏差统计结果图。从本发明方法与UC方法共同收敛部分(3007组数据)的收敛时间和平面位置偏差统计结果可以得到：本发明方法和UC方法的平均收敛时间分别为12.3min和13.6min，较UC方法缩短了9.6％，且分别有50％的收敛时间小于11.0min和15.5min；本发明方法和UC方法估计得到的接收机平面位置的平均偏差分别为3.7cm和4.3cm，且分别有75.4％和65.0％的平面位置偏差位于(0，5)cm内。

步骤(3)：选取本发明方法与UofC方法公共滤波收敛部分，使用本发明方法与UofC方法对收敛时间和平面位置偏差方面进行分析，得到收敛时间统计结果和平面位置偏差统计结果，参见图5和图6，图5为使用本发明方法和UofC方法对公共滤波收敛数据处理后的收敛时间统计结果图；图6为使用本发明方法和UofC方法对公共滤波收敛数据处理后的平面位置偏差统计结果图。从本发明方法与UofC方法共同收敛部分(4090组数据)的收敛时间和平面位置偏差统计结果可以得到：本发明方法和UofC方法的平均收敛时间分别为12.1min和15.7min，较UofC方法缩短了22.9％，且分别有50％的收敛时间小于11.5min和13.0min；本发明方法和UofC方法估计得到的接收机平面位置的平均偏差分别为3.7cm和3.9cm，且分别有73.9％和71.5％的平面位置偏差位于(0，5)cm内。

步骤(4)：选取本发明方法与UD方法公共滤波收敛部分，使用本发明方法与UD方法对收敛时间和平面位置偏差方面进行分析，得到收敛时间统计结果和平面位置偏差统计结果，参见图7和图8，图7为使用本发明方法和UD方法对公共滤波收敛数据处理后的收敛时间统计结果图；图8为使用本发明方法和UD方法对公共滤波收敛数据处理后的平面位置偏差统计结果图。从本发明方法与UD方法共同收敛部分(1180组数据)的收敛时间和平面位置偏差统计结果可以得到：本发明方法和UD方法的平均收敛时间分别为11.9min和18.4min，较UD方法缩短了35.3％，且分别有50％的收敛时间小于11.0min和19.0min；本发明方法和UD方法估计得到的接收机平面位置的平均偏差分别为3.7cm和5.2cm，且分别有74.7％和48.2％的平面位置偏差位于(0，5)cm内。

通过本应用实施例的数据处理，得到：使用本发明的一种快速收敛的精密单点定位方法能够快速、准确、高比例的估计接收机位置信息，相对于目前存在的UC方法、UofC方法和UD方法在定位精度和收敛速度方面均有明显提高。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同技术特征的含义和范围内的所有变化囊括在本发明内。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (3)

1.一种快速收敛的精密单点定位方法,其特征在于,包括如下步骤:

(S1)、对全球导航卫星的观测数据进行预处理：

采用M-W组合法和电离层残差法对原始载波相位观测值中存在的周跳进行探测，并修复；

对天线相位中心偏差、地球自转、相对论效应、天线相位缠绕误差源的影响进行修正；

对精密卫星轨道和卫星钟差数据进行拉格朗日多项式内插，获取观测时刻的卫星轨道和卫星钟差；

(S2)、顾及电离层延迟对不同频率观测值的影响关系，建立卫星的原始载波相位观测值和的观测方程，如公式(5)、公式(6)所示：

${\mathrm{\Φ}}_1^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_1^s}---\left(5\right)$

${\mathrm{\Φ}}_2^s=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}-f_1^2/f_2^2\·I_1^s+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_2^s}---\left(6\right)$

其中，s表示卫星，和表示载波相位观测值，ρ表示卫星与接收机之间的几何距离，c表示真空中光的速度，表示接收机钟差，表示卫星钟差，分别表示对流层干、湿迟延映射函数；δ_zhd、δ_zwd分别表示对流层天顶方向干延迟和湿延迟，为载波相位观测值中站星视线方向的电离层延迟，分别为原始载波相位观测值和的模糊度，f₁、f₂为载波相位的频率，λ₁、λ₂表示载波相位的波长，表示载波相位观测值的观测噪声；

(S3)、根据原始码伪距观测值和原始载波相位观测值构建码相位半合组合观测值的观测方程，如公式(7)所示：

$L_{P_2{\mathrm{\Φ}}_2}^s=(P_2^s+{\mathrm{\Φ}}_2^s)/2=\mathrm{\ρ}+c\·\mathrm{\δ}\stackrel{\‾}{t_r}-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zpd}+M_w^s\·{\mathrm{\δ}}_{zpw}+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s/2+{\mathrm{\ϵ}}_{L_{P_2{\mathrm{\Φ}}_2}^s}---\left(7\right)$

(S4)、根据原始载波相位观测值和构建几何无关组合观测值的观测方程，如公式(8)所示：

${\mathrm{\Φ}}_{GF}^s={\mathrm{\Φ}}_1^s-{\mathrm{\Φ}}_2^s=(f_1^2/f_2^2-1)\·I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s-{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_{GF}^s}---\left(8\right)$

(S5)、联合原始载波相位观测值和的观测方程、码相位半合组合观测值的观测方程和几何无关组合观测值的观测方程构建改进的精密单点定位的观测方程，如公式(9)所示，

$\left\{\begin{array}{c}{\mathrm{\Φ}}_1^s=\mathrm{\ρ}+c\·\mathrm{\δ}{\stackrel{\‾}{t}}_r-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zqd}+M_w^s\·{\mathrm{\δ}}_{zqw}-I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_2^s}\\ {\mathrm{\Φ}}_2^s=\mathrm{\ρ}+c\·\mathrm{\δ}{\stackrel{\‾}{t}}_r-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zqd}+M_w^s\·{\mathrm{\δ}}_{zqw}-f_1^2/f_2^2\·I_1^s+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_2^s}\\ L_{P_2{\mathrm{\Φ}}_2}^s=(P_2^s+{\mathrm{\Φ}}_2^s)/2=\mathrm{\ρ}+c\·\mathrm{\δ}{\stackrel{\‾}{t}}_r-c\·\mathrm{\δ}{\stackrel{\‾}{t}}^s+M_d^s\·{\mathrm{\δ}}_{zqd}+M_w^s\·{\mathrm{\δ}}_{zqw}+{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s/2+{\mathrm{\ϵ}}_{L_{P_2{\mathrm{\Φ}}_2}^s}\\ {\mathrm{\Φ}}_{GF}^s={\mathrm{\Φ}}_1^s-{\mathrm{\Φ}}_2^s=(f_1^2/f_2^2-1)\·I_1^s+{\mathrm{\λ}}_1\·{\stackrel{\‾}{N}}_1^s-{\mathrm{\λ}}_2\·{\stackrel{\‾}{N}}_2^s+{\mathrm{\ϵ}}_{{\mathrm{\Φ}}_{GF}^s}\end{array}\right.---\left(9\right)$

(S6)、采用卡尔曼滤波对改进的精密单点定位的观测方程中的参数进行估计，得到接收机位置信息和对流层延迟、电离层延迟信息。

2.根据权利要求1所述的一种快速收敛的精密单点定位方法,其特征在于:对步骤(S5)的改进的精密单点定位的观测方程进行线性化处理，得到公式(10)：

y(k)＝A·X(k)+ε_y，ε_y～N(0,Ω_y)(10)

其中，y(k)为观测向量，其维数为4j，j为历元k同步观测到的卫星数量；X(k)为待估参数向量，维数为(5+3j)，A为设计矩阵，ε_y为测量噪声向量，Ω_y为测量噪声向量ε_y的协方差阵；y(k)、X(k)和A表达式分别为

$y\left(k\right)={\left[\begin{array}{cccccccccccc}{\mathrm{\Φ}}_1^1& ...& {\mathrm{\Φ}}_1^j& {\mathrm{\Φ}}_2^1& ...& {\mathrm{\Φ}}_2^j& L_{P_2{\mathrm{\Φ}}_2}^1& ...& L_{P_2{\mathrm{\Φ}}_2}^j& {\mathrm{\Φ}}_{GF}^1& ...& {\mathrm{\Φ}}_{GF}^j\end{array}\right]}^T$

$X\left(k\right)={\left[\begin{array}{cccccccccccccc}X& Y& Z& c\·\mathrm{\δ}\stackrel{\‾}{t_r}& {\mathrm{\δ}}_{zpw}& I_1^1& ...& I_1^j& {\stackrel{\‾}{N}}_1^1& ...& {\stackrel{\‾}{N}}_1^j& {\stackrel{\‾}{N}}_2^1& ...& {\stackrel{\‾}{N}}_2^j\end{array}\right]}^T$

$A=\[\begin{array}{ccc}{e}_4\⊗B& C\⊗I_j& D\⊗I_j\end{array}\]$

其中，X、Y、Z为接收机三维位置，表示克罗内克积，e₄为各元素均为1的4维列向量，I_j为j维单位矩阵，B、C和D子矩阵的含义分别为:

$B=\left[\begin{array}{ccccc}{l}^1& m^1& n^1& M_w^1& 1\\ \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}& \begin{array}{c}\·\\ \·\\ \·\end{array}\\ l^j& m^j& n^j& M_w^j& 1\end{array}\right],C=\left[\begin{array}{c}-1\\ -f_1^2/f_2^2\\ 0\\ f_1^2/f_2^2-1\end{array}\right],D=\left[\begin{array}{cc}{\mathrm{\λ}}_1& 0\\ 0& {\mathrm{\λ}}_2\\ 0& {\mathrm{\λ}}_2/2\\ {\mathrm{\λ}}_1& -{\mathrm{\λ}}_2\end{array}\right],$

其中，l、m、n分别为接收机到卫星的方向余弦。

3.根据权利要求1所述的一种快速收敛的精密单点定位方法,其特征在于:步骤(S5)得到的改进的精密单点定位的观测方程第k历元对应的协方差阵Ω_y(k)，如公式(11)所示：

${\mathrm{\Ω}}_y\left(k\right)=F\⊗I_j---\left(11\right)$

其中， $F=\left[\begin{array}{cccc}{\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2& 0& 0& {\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2\\ 0& {\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2& \frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}& -{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2\\ 0& \frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}& \frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2+{\mathrm{\σ}}_{P_2}^2}{4}& -\frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}\\ {\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2& -{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2& -\frac{{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2}{2}& {\mathrm{\σ}}_{{\mathrm{\Φ}}_1}^2+{\mathrm{\σ}}_{{\mathrm{\Φ}}_2}^2\end{array}\right],$ 表示克罗内克积，I_j为j维单位矩阵，分别为卫星原始载波相位观测值Φ₁、Φ₂，原始码伪距观测值P₂的测量噪声的方差。