CN101403790B - 单频gps接收机的精密单点定位方法 - Google Patents

Abstract

本发明涉及一种单频GPS接收机的精密单点定位方法，包括GPS数据的读取和预处理步骤、误差改正解算步骤和结果输出步骤，所述预处理步骤包括进行周跳的探测和修复，所述周跳的探测和修复步骤包括：1)将m个无周跳的载波相位观测值φ_i带入公式，进行多项式拟合；2)用求得的多项式系数来外推下一个历元的载波相位观测值

，并与实际观测值φ_i进行比较；3)利用第i和i-1历元载波相位观测值，获得载波相位观测值的变化率。本发明综合利用了多项式拟合法和载波相位变化率法来探测周跳，结合粗差修改实现单频GPS接收机的精密单点定位，有效实现了GPS系统的高精度定位和导航。

Description

单频GPS接收机的精密单点定位方法

(一)技术领域

本发明涉及一种单频GPS接收机的精密单点定位方法。

(二)背景技术

GPS全球定位系统是英文“Navigation Satellite Timing And Ranging/GlobalPositioning System”的头字母缩写词NAVSTAR/GPS的简称。它是二十世纪七十年代美国国防部开始研制的新一代卫星导航系统，是基于卫星的无线电导航定位系统，具有全能性(陆地、海洋、航空和航天)，全球性、全天候、连续性和实时性的导航、定位和授时功能，能为各类用户提供精密的三维定位、测速和时间服务。目前，GPS作为新一代的卫星导航定位系统已经在军事、交通运输、测绘、高精度时间比对、以及资源调查等领域得到了广泛的应用；它在动态导航、天气预报等方面的应用与人们的日常生活息息相关；同时在航海导航、飞机导航、时间转换、救助、制图、导弹导航等复杂工作中，GPS也是不可缺少的。据专家预测，今后几年内GPS在通信、大气探测、精细农业以及环境保护等领域中也将会得到广泛的应用。

GPS系统主要由空间部分、地面控制部分、用户部分组成。卫星播放的信号包括L1、C\A码、P码、L2、P码。对于单频接收机用户，能获得的信号只有L1和C\A码。GPS定位模式主要有非差和差分模式。非差模式又分为伪距单点定位和相位精密单点定位；差分模式分为单差、双差和三差等模式。

伪距单点定位是指利用伪距码和广播星历，采用距离交会法解算接收机天线所在点的三维坐标。伪距单点定位的数据采集和数据处理简单，定位速度快，用户在任一时刻只需使用一台GPS接收机就能获得测点在WGS-84坐标系中的三维坐标。但这种单点定位的精度只能达到十米级，只能满足一些低精度导航定位领域的需求。

GPS导航定位技术通过几十年的变化发展，提出了许多获得厘米级精度的定位方法。目前广泛使用的是差分GPS定位方法，通过组成双差观测值消除接收机钟差、卫星钟差等公共误差，削弱对流层延迟、电离层延迟等相关性强的误差影响，达到提高精度的目的。差分GPS无需考虑复杂的误差模型，解算模型简单、待估参数少、定位精度高，同时利用了双差模糊度的整数特征，因而得到广泛使用。其不足之处在于：1)至少需要使用两台GPS接收机进行同步观测，一个参考站和一个移动站；2)移动站和参考站之间必须存在较强的空间关联；3)参考站坐标误差会直接传给移动站坐标，用户在几个精度不同的参考站间转换，定位结果就会不一致。

针对差分GPS的不足和传统GPS单点定位的低精度，JPL的Zumbeger等人于1997年首先提出了非差精密单点定位方法，即利用高精度的GPS卫星星历和卫星钟差以及双频载波相位观测值，采用非差模型进行精密单点定位。精密单点定位技术克服了传统差分GPS的缺点，根据IGS和其它机构提供的精密星历和卫星钟差消除由轨道和钟差引起的误差；使用双频GPS观测值消除电离层延迟；使用精确非差载波相位观测值减弱其它误差项，从而得到高精度的绝对定位结果，提高了GPS精密定位的操作机动性，减少了投资成本。

近几年来，单频GPS精密定位技术吸引了很多研究者进行研究。其主要原因是硬件的廉价和后处理方面的优势。单频精密单点定位(Single-frequencyPrecise Point Positioning)是最近几年发展起来的一项GPS定位新技术，它利用IGS或者其他机构发布的卫星精密星历及由一定算法确定的精密卫星钟差，以单台单频GPS接收机采集的相位观测数据和伪距观测数据作为主要观测值来进行单点定位计算，在全球范围内其精度可达亚米级。由于它可利用单台单频接收机在全球范围内进行静态或动态独立作业，并可以直接得到高精度的ITRF框架坐标，因此，作业人员只需要携带单频GPS接收机即可方便地解决定位问题；除此之外在移动电话中配备GPS接收芯片，用户可以很方便地得到准确的定位结果；当然在高精度动态导航定位及低轨卫星的定轨等方面都具有很好的应用前景。

双频非差精密单点定位研究已经相当成熟，不过，考虑到双频接收机的价格昂贵，以及某些应用未必需要到达厘米或毫米级的精度。所以，近几年来，单频精密单点定位技术吸引了很多国内外研究员的研究。其中包括挪威学者Ola

荷兰学者Anh Quan Le，加拿大学者Tomas Beran等，他们都发表了关于单频精密单点定位的文章并且详细阐述了其中的技术细节。由于单频接收机只能接收L1载波相位观测数据和C/A码伪距。单频精密单点定位相比双频精密单点定位有如下难点：

1)单点定位中部分误差无法采用求差的方式消除，必须使用各种模型和参数估计方法进行改正：

a.对流层湿分量的不规则性使得模型的改正精度比较低，需要使用参数估计的方法来确定；

b.无法采用双差观测值来消除电离层的影响，必须使用电离层模型来减弱它的影响，或者利用伪距和相位组合的半和改正；

2)观测值中噪声很较大，不能通过组合值的方法探测，必然加大了周跳探测的难度，特别是对于含有小周跳的观测值；

3)由于单站数据无法采用双差观测值来消除各种误差，模型改正无法完全消除各项误差，而且非差相位观测值的模糊度本身就不一定具有整数特性，所以其模糊度解算非常困难。

(三)发明内容

本发明所要解决的技术问题在于提供一种在GPS数据预处理步骤中正确探测与修复周跳，结合有效的误差改正，从而提高定位、定轨或导航精度的单频GPS接收机的精密单点定位方法。

所述的单频GPS接收机的精密单点定位方法，包括GPS数据的读取和预处理步骤、误差改正解算步骤和结果输出步骤，所述预处理步骤包括进行周跳的探测和修复，所述周跳的探测和修复步骤包括：

1)将m个无周跳的载波相位观测值φ_i带入公式(1)，进行多项式拟合；

用最小二乘法求得式中的多项式系数，并根据拟合后的残差V_i计算出误差：

$\mathrm{\σ}\left(j\right)=\sqrt{\frac{V_i{V}_i}{m-(n+1)}},$ 其中j表示第j次拟合                (2)

2)用求得的多项式系数来外推下一个历元的载波相位观测值

并与实际观测值φ_i进行比较：

其中k为门阀系数                   (3)

若公式(3)成立，则认为该观测值没有周跳；

加入上述无周跳的实际观测值后继续上述过程进行多项式拟合；

若公式(4)成立，则认为实际观测值发生了周跳，此时应采用外推的整周计数去取代有周跳的实际观测值中的整周计数，但不足一周的部分F_r(φ)仍保持不变；继续上述过程，直到最后一个观测值为止；

3)利用第i和i-1历元载波相位观测值，获得载波相位观测值的变化率：

                      φ′_i＝(φ_i-φ_i-1)/Δt                (5)

又由：

                      φ′_i+1＝φ′_i                        (6)

基于第i历元的观测值，并利用上式(6)，可预报第i+1历元的观测值：

若载波相位不出现周跳，它与第i+1历元的观测值的差值应在一定范围之内：

如果公式(8)不成立，则认为有周跳发生并进行修复，精度可在±2周之内。

所述的步骤3)中运用载波相位变化率探测并修复周跳，公式(6)的依据在于：由于频率间偏差和电离层延迟的变化都非常缓慢，载波相位观测模糊度不随时间变化，其测值噪声和多路径效应变化幅度又很小，所以Φ′_i基本反映的是伪距的变化率(包含接收机钟差和卫星钟差的变化率)。同时，卫星相对于静止测站的速度、接收机钟速和卫星钟速都很稳定，即伪距变化率在短时间内变化很小。

优选地，所述的误差改正解算步骤包括：

4)把L1和C/A码进行组合得到一个新的观测值La，并建立以La为观测值的观测方程，若测站的近似坐标为(x_r ⁰、y_r ⁰、z_r ⁰)，将观测方程在(x_r ⁰、y_r ⁰、z_r ⁰)处用泰勒级数展开得到线性化的观测方程。将误差方程与同一历元中GPS卫星i的观测方程相减，即以GPS卫星i作为参考卫星，得到：

$V_{\mathrm{ij}}=-(l_j-l_i)V_x-(m_j-m_i)V_y-(n_j-n_i)V_z-{\mathrm{\λ}}_1(N_a^j-N_a^i)-C({\mathrm{dt}}^j-{\mathrm{dt}}^i)+(R^j-R^i)---\left(10\right)$

从而消除了电离层一阶项的影响和接收机钟差的影响；

5)卫星质量中心与卫星发射天线的相位中心的不重合产生天线相位中心偏差：该偏差向量在常用的星固坐标系中的三个坐标分量为：星固系的Z轴指向地心，其单位向量为e_z；星固系的Y轴是卫星方向与太阳方向至卫星方向的叉乘，其单位向量为e_y；星固系的X轴与另外两个轴组成右手系，其单位向量为e_x。因此卫星坐标系从质量中心到相位中心的改正值为：

$\mathrm{\Δ}{R}_{\mathrm{sat}}=(e_x,e_y,e_z)\left(\begin{array}{c}{\mathrm{\ΔX}}_{\mathrm{sat}}\\ {\mathrm{\ΔY}}_{\mathrm{sat}}\\ {\mathrm{\ΔZ}}_{\mathrm{sat}}\end{array}\right)---\left(11\right)$

6)利用伪距做初始解算，获得接收机的近似位置；

7)采用Hopfield模型改正天顶方向对流层延迟：

${\mathrm{\ΔP}}_d=155.2\×{10}^{-7}\×\frac{P_s}{T_s}\×(h_d-h_s)---\left(12\right)$

${\mathrm{\ΔP}}_w=155.2\×{10}^{-7}\×\frac{4810}{{T_s}^2}\×e_s\×(h_w-h_s)---\left(13\right)$

                     h_d＝40136+148.72×(T_s-273.16)              (14)

                               h_w＝11000                        (15)

             e_s＝RH·exp(-37.2465+0.21366T-0.0002568908T²)      (16)

其中，T_s，P_s，e_s，h_s分别为测站的温度、气压、湿度、正高，RH为相对湿度，T为干温，为开尔文温度；

8)卫星天线在进出地影时，太阳能面板会慢慢地旋转，导致卫星天线重新定向，改变测站和卫星之间的几何关系，产生相位缠绕，改正方法为：

$\mathrm{\ζ}=\hat{k}\·({\stackrel{\⊥}{D}}^{\′}\·\stackrel{\⊥}{D})---\left(18\right)$

${\stackrel{\⊥}{D}}^{\′}={\hat{x}}^{\′}-\hat{k}(\hat{k}\·{\hat{x}}^{\′})-\hat{k}\×{\hat{y}}^{\′}---\left(19\right)$

$\stackrel{\⊥}{D}=\hat{x}-\hat{k}(\hat{k}\·\hat{x})-\hat{k}\×\hat{y}---\left(20\right)$

其中，是卫星到接收机的单位向量，

是由星固坐标系下的单位向量

计算得到的卫星有效的偶极向量，

是由地方坐标系下的单位向量

计算得到的接收机有效的偶极向量；

9)利用拉格朗日插值算法对精密星历和精密钟差进行任意时刻的插值：设在n+1个时间为t₁，t₂，...t_n+1插值节点上卫星的坐标分别是x(t₁)，x(t₂)……x(t_n+1)。那么在任意时刻t卫星的坐标用下述公式表示：

$x\left(t\right)=\underset{i=1}{\overset{n+1}{\∑}}[x\left(t_i\right)\·\underset{j\≠i}{\overset{n+1}{\underset{j=1}{\mathrm{\Π}}}}\frac{t-t_j}{t_i-t_j}]---\left(21\right)$

利用公式在x、y和z方向上以及根据精密钟差文件分别进行插值处理，即可得到任意时刻卫星的位置和钟差值；

10)万有引力作用下使得地球表面产生固体潮现象，固体潮对测站位置的影响公式如下：

其中：r_E为地球的半径；X_j为摄动天体在地心参考框架中的坐标向量；X_p为测站在地心参考框架中的坐标向量；GM_i为摄动天体(j＝2为月球，j＝3为太阳)的引力参数；GM_E为地球引力参数；h₂、l₂为Love和Shida数(h₂＝0.6090与l₂＝0.0852)；φ，λ为测站纬度和经度；θ为格林尼治平恒星时；

11)运用卡尔曼滤波法，利用受到噪声污染的观测值，对对流层参数、模糊度参数做最优化的估计；

所述的结果输出步骤包括：

12)输出地固系下的(X，Y，Z)或者是大地坐标(B，L，H)，以及在以平均坐标(或者已知真实坐标)为原点的站心地平坐标下的(N，E，U)等计算结果。

优选地，所述的步骤1)中，公式(1)中n＝4，m＝14。由于卫地距对时间的四阶导数或五阶导数一般已趋近于零了，其变化规律是随机的，无法再用多项式拟合，因此阶数n取4阶即可；当拟合窗宽度越大，虽然外推的值越准确，同时拟合后的中误差会很小，导致公式(3)很难满足，而且计算量会加大。但是当拟合窗宽度越小，外推的值会越粗糙，同样也导致公式(3)很难满足。通过取不同的值进行试验，在周跳探测时，取m＝14最为优选。

优选地，所述的步骤2)中，公式(3)中k＝5。式(3)中的k值是根据情况而设定的。当k较大时，表示只有当观测值偏离外推值很大的时候才认为它是异常的；当k较小时，表示当观测值偏离外推值较小的时候就认为它是异常值了。一般取k在2到9之间的数。经过发明人多次取值实验，k＝5最为优选。

优选地，所述的步骤4)中利用半合改正来消除电离层影响，包括如下具

(4.1)把L1和C/A码进行组合得到一个新的观测值：

其模糊度为则以La为观测值的观测方程为：

$L_a^j={\mathrm{\ρ}}^j-{\mathrm{\λ}}_1{N}_a^j-C({\mathrm{dt}}^j-\mathrm{dT})+{\mathrm{dtrop}}^j+{\mathrm{dmult}}_a^j+{\mathrm{drel}}^j+{\mathrm{\ϵ}}_a^j---\left(23\right)$

(4.2)将公式(23)在(x_r ⁰、y_r ⁰、z_r ⁰)处用泰勒级数展开后，得到线性化的观测方程：

$L_a^j={\mathrm{\ρ}}_0^j-\frac{x_s^j-x_r^0}{{\mathrm{\ρ}}_0^j}{V}_x-\frac{y_s^j-y_r^0}{{\mathrm{\ρ}}_0^j}{V}_y-\frac{z_s^j-z_r^0}{{\mathrm{\ρ}}_0^j}{V}_z$

$-{\mathrm{\λ}}_1{N}_a^j-C({\mathrm{dt}}^j-\mathrm{dT})+{\mathrm{dtrop}}^j+{\mathrm{dmult}}_{\mathrm{cA}}^j+{\mathrm{drel}}^j+{\mathrm{\ϵ}}_{\mathrm{cA}}^j---\left(24\right)$

其中，令： $\frac{x_s^j-x_r^0}{{\mathrm{\ρ}}_0^j}=l_j,$ $\frac{y_s^j-y_r^0}{{\mathrm{\ρ}}_0^j}=m_j,$ $\frac{z_s^j-z_r^0}{{\mathrm{\ρ}}_0^j}=n_j$

(4.3)误差方程表示为：

$V_j=-l_j{V}_x-m_j{V}_y-n_j{V}_z-{\mathrm{\λ}}_1{N}_a^j-C({\mathrm{dt}}^j-\mathrm{dT})+R^j---\left(25\right)$

其中： $R^j={\mathrm{\ρ}}_0^j-L_a^j+d_{{\mathrm{trop}}^j}+{\mathrm{drel}}^j$

对于同一历元，对GPS卫星i的观测方程为：

$V_i=-l_i{V}_x-m_i{V}_y-n_i{V}_z-{\mathrm{\λ}}_1{N}_a^i-C({\mathrm{dt}}^i-\mathrm{dT})+R^i---\left(26\right).$

优选地，所述的步骤5)中天线相位中心偏差改正，包括如下具体步骤：

星固系的Z轴指向地心，其单位向量e_z为：

$e_z=\frac{-r_{\mathrm{sat}}}{{|r}_{\mathrm{sat}}|}---\left(27\right)$

式中的r_sat为卫星质量中心的坐标；

星固系的Y轴是卫星方向与太阳方向至卫星方向的叉乘，其单位向量e_y为：

$e_y=\frac{r_{\mathrm{sat}}\×(r_{\mathrm{sat}}-r_{\mathrm{sun}})}{|r_{\mathrm{sat}}\×(r_{\mathrm{sat}}-r_{\mathrm{sun}})|}---\left(28\right)$

式中的r_sun为太阳坐标；

星固系的X轴与另外两个轴组成右手系，其单位向量e_x为：

星固系的X轴与另外两个轴组成右手系，其单位向量e_x为：

$e_x=\frac{e_y\×e_z}{|e_y\×e_z|}---\left(29\right).$

优选地，所述的步骤11)中运用卡尔曼滤波进行参数估计，包括如下具体步骤：

(11.1)假设线性离散系统的状态方程和观测方程为：

            X_k+1＝Φ_k+1，kX_k+Г_kω_k                        (30)

            Z_k+1＝H_k+1X_k+1+v_k+1                            (31)

其中，X_k为n维状态向量；Φ_k+1为n×n维的一步状态转移矩阵；Γ_k为n×p维动态噪声驱动阵，ω_k为P维系统动态噪声向量，Z_k+1为m维观测向量，H_k+1为m×n观测矩阵，v_k+1为m维观测噪声向量；

(11.2)首先根据前一次滤波值

(或初值)计算预测值：

$\hat{X}(k/k-1)=\mathrm{\Φ}(k,k-1)\hat{X}(k-1/k-1)---\left(32\right)$

根据前一次得到的滤波误差方差阵P(k-1/k-1)(或初值)及系统噪声的方差阵Q_k计算预测误差方差阵：

P(k/k-1)＝Q(k，k-1)P(k-1/k-1)Q^T(k/k-1)-Г(k-1)Q(k-1)Г^T(k-1)    (33)

(11.3)计算卡尔曼滤波增益，公式如下：

K(k)＝P(k-1/k-1)H^T(k)[H(k)P(k/k-1)H^T(k)+R(k)]^T                  (34)

根据新的观测值z(k)得到：

$v\left(k\right)=z\left(k\right)-H\left(k\right)\hat{X}(k/k-1)---\left(35\right)$

(11.4)计算滤波估计，公式如下：

$\hat{X}(k/k)=\hat{X}(k-1/k-1)+K\left(k\right)v\left(k\right)---\left(36\right)$

计算滤波误差方差阵，公式如下：

                P(k/k)＝[I-K(k)H(k)]P(k/k-1)                    (37)

(11.5)将滤波估计存入计算机，下一时刻得到新的观测值，重复上述计算过程，滤波以不断“预测-修正”的递推方式进行计算，对预测值不断进行修正。

本发明综合利用了多项式拟合法和载波相位变化率法来探测周跳，结合粗差修改实现单频GPS接收机的精密单点定位。多项式拟合法是以星间单差作为处理对象，可探测出组合观测值1周以上的周跳。载波相位变化率法可判断是哪一颗卫星发生了周跳。二者相互补充，可以正确探测与修复周跳。同时，误差改正解算步骤中采用半和改正法，由于消除了电离层的影响，其定位精度在N、E、U方向都能达到2分米左右。上述技术方案有效实现GPS系统的高精度定位和导航。

(四)附图说明

图1为单差观测值求四阶差分的时间序列图。

图2为每次拟合的中误差图。

图3为外推值与实际观测值之差图。

图4为外推值跟实际观测值之差与中误差的比值图。

图5为在23号卫星的L1上的第25个历元模拟加上3个周跳图。

图6为载波相位变化率法探测23号卫星(选择了前100个历元)图。

图7为载波相位变化率法探测13号卫星(选择了前100个历元)图。

图8为格网模型改正解算结果在N、E、U方向统计图(bjfs0150.07O)。

图9为格网模型改正解算结果在N、E、U方向统计图(wuhn0150.07O)。

图10为半合改正解算结果在N、E、U方向统计图(bjfs0150.07O)。

图11为半合改正解算结果在N、E、U方向统计图(wuhn0150.07O)。

图12为单频GPS精密单点定位方法的流程图。

(五)具体实施方式

下面通过实施例对本发明作优选地具体的说明，但本发明的保护范围并不限于此。

实施例1：

参照图12，一种单频GPS接收机的精密单点定位方法，包括GPS数据的读取和预处理步骤、误差改正解算步骤和结果输出步骤，所述预处理步骤包括进行周跳的探测和修复，所述周跳的探测和修复步骤包括：

1)将m个无周跳的载波相位观测值φ_i带入公式(1)，进行多项式拟合；

用最小二乘法求得式中的多项式系数，并根据拟合后的残差V_i计算出误差：

$\mathrm{\σ}\left(j\right)=\sqrt{\frac{V_i{V}_i}{m-(n+1)}},$ 其中j表示第j次拟合                (2)

2)用求得的多项式系数来外推下一个历元的载波相位观测值

并与实际观测值φ_i进行比较：

其中k为门阀系数                   (3)

若公式(3)成立，则认为该观测值没有周跳；

加入上述无周跳的实际观测值后继续上述过程进行多项式拟合；

若公式(4)成立，则认为实际观测值发生了周跳，此时应采用外推的整周计数去取代有周跳的实际观测值中的整周计数，但不足一周的部分F_r(φ)仍保持不变；继续上述过程，直到最后一个观测值为止；

3)利用第i和i-1历元载波相位观测值，获得载波相位观测值的变化率：

                     φ′_i＝(φ_i-φ_i-1)/Δt                (5)

又由：

                     φ′_i+1＝φ′_i                        (6)

基于第i历元的观测值，并利用上式(6)，可预报第i+1历元的观测值：

若载波相位不出现周跳，它与第i+1历元的观测值的差值应在一定范围之内：

如果公式(8)不成立，则认为有周跳发生并进行修复，精度可在±2周之内。

所述的误差改正解算步骤包括：

4)把L1和C/A码进行组合得到一个新的观测值La，并建立以La为观测值的观测方程，若测站的近似坐标为(x_r ⁰、y_r ⁰、z_r ⁰)，将观测方程在(x_r ⁰、y_r ⁰、z_r ⁰)处用泰勒级数展开得到线性化的观测方程。将误差方程与同一历元中GPS卫星i的观测方程相减，即以GPS卫星i作为参考卫星，得到：

$V_{\mathrm{ij}}=-(l_j-l_i)V_x-(m_j-m_i)V_y-(n_j-n_i)V_z-{\mathrm{\λ}}_1(N_a^j-N_a^i)-C({\mathrm{dt}}^j-{\mathrm{dt}}^i)+(R^j-R^i)---\left(10\right)$

从而消除了电离层一阶项的影响和接收机钟差的影响；

5)卫星质量中心与卫星发射天线的相位中心的不重合产生天线相位中心偏差：该偏差向量在常用的星固坐标系中的三个坐标分量为：星固系的Z轴指向地心，其单位向量为e_z；星固系的Y轴是卫星方向与太阳方向至卫星方向的叉乘，其单位向量为e_y；星固系的X轴与另外两个轴组成右手系，其单位向量为e_x。因此卫星坐标系从质量中心到相位中心的改正值为：

$\mathrm{\Δ}{R}_{\mathrm{sat}}=(e_x,e_y,e_z)\left(\begin{array}{c}{\mathrm{\ΔX}}_{\mathrm{sat}}\\ {\mathrm{\ΔY}}_{\mathrm{sat}}\\ {\mathrm{\ΔZ}}_{\mathrm{sat}}\end{array}\right)---\left(11\right)$

6)利用伪距做初始解算，获得接收机的近似位置；

7)采用Hopfield模型改正天顶方向对流层延迟：

${\mathrm{\ΔP}}_d=155.2\×{10}^{-7}\×\frac{P_s}{T_s}\×(h_d-h_s)---\left(12\right)$

${\mathrm{\ΔP}}_w=155.2\×{10}^{-7}\×\frac{4810}{{T_s}^2}\×e_s\×(h_w-h_s)---\left(13\right)$

                      h_d＝40136+148.72×(T_s-273.16)            (14)

                                h_w＝11000                      (15)

              e_s＝RH·exp(-37.2465+0.21366T-0.0002568908T²)    (16)

其中，T_s，P_s，e_s，h_s分别为测站的温度、气压、湿度、正高，RH为相对湿度，T为干温，为开尔文温度；

8)卫星天线在进出地影时，太阳能面板会慢慢地旋转，导致卫星天线重新定向，改变测站和卫星之间的几何关系，产生相位缠绕，改正方法为：

$\mathrm{\ζ}=\hat{k}\·({\stackrel{\⊥}{D}}^{\′}\·\stackrel{\⊥}{D})---\left(18\right)$

${\stackrel{\⊥}{D}}^{\′}={\hat{x}}^{\′}-\hat{k}(\hat{k}\·{\hat{x}}^{\′})-\hat{k}\×{\hat{y}}^{\′}---\left(19\right)$

$\stackrel{\⊥}{D}=\hat{x}-\hat{k}(\hat{k}\·\hat{x})-\hat{k}\×\hat{y}---\left(20\right)$

其中，

是卫星到接收机的单位向量，

是由星固坐标系下的单位向量

计算得到的卫星有效的偶极向量，是由地方坐标系下的单位向量

计算得到的接收机有效的偶极向量；

9)利用拉格朗日插值算法对精密星历和精密钟差进行任意时刻的插值：设在n+1个时间为t₁，t₂，...t_n+1插值节点上卫星的坐标分别是x(t₁)，x(t₂)……x(t_n+1)。那么在任意时刻t卫星的坐标用下述公式表示：

$x\left(t\right)=\underset{i=1}{\overset{n+1}{\∑}}[x\left(t_i\right)\·\underset{j\≠i}{\overset{n+1}{\underset{j=1}{\mathrm{\Π}}}}\frac{t-t_j}{t_i-t_j}]---\left(21\right)$

利用公式在x、y和z方向上以及根据精密钟差文件分别进行插值处理，即可得到任意时刻卫星的位置和钟差值；

10)万有引力作用下使得地球表面产生固体潮现象，固体潮对测站位置的影响公式如下：

其中：r_E为地球的半径；X_j为摄动天体在地心参考框架中的坐标向量；X_p为测站在地心参考框架中的坐标向量；GM_i为摄动天体(j＝2为月球，j＝3为太阳)的引力参数；GM_E为地球引力参数；h₂、l₂为Love和Shida数(h₂＝0.6090与l₂＝0.0852)；φ，λ为测站纬度和经度；θ为格林尼治平恒星时；

11)运用卡尔曼滤波法，利用受到噪声污染的观测值，对对流层参数、模糊度参数做最优化的估计；

所述的结果输出步骤包括：

12)输出地固系下的(X，Y，Z)或者是大地坐标(B，L，H)，以及在以平均坐标(或者已知真实坐标)为原点的站心地平坐标下的(N，E，U)等计算结果。

所述的步骤1)中，公式(1)中n＝4，m＝14。

所述的步骤2)中，公式(3)中k＝5。

所述的步骤4)中利用半合改正来消除电离层影响，包括如下具体步骤：

(4.1)把L1和C/A码进行组合得到一个新的观测值：

其模糊度为

则以La为观测值的观测方程为：

$L_a^j={\mathrm{\ρ}}^j-{\mathrm{\λ}}_1{N}_a^j-C({\mathrm{dt}}^j-\mathrm{dT})+{\mathrm{dtrop}}^j+{\mathrm{dmult}}_a^j+{\mathrm{drel}}^j+{\mathrm{\ϵ}}_a^j---\left(23\right)$

(4.2)将公式(23)在(x_r ⁰、y_r ⁰、z_r ⁰)处用泰勒级数展开后，得到线性化的观测方程：

$L_a^j={\mathrm{\ρ}}_0^j-\frac{x_s^j-x_r^0}{{\mathrm{\ρ}}_0^j}{V}_x-\frac{y_s^j-y_r^0}{{\mathrm{\ρ}}_0^j}{V}_y-\frac{z_s^j-z_r^0}{{\mathrm{\ρ}}_0^j}{V}_z$

$-{\mathrm{\λ}}_1{N}_a^j-C({\mathrm{dt}}^j-\mathrm{dT})+{\mathrm{dtrop}}^j+{\mathrm{dmult}}_{\mathrm{cA}}^j+{\mathrm{drel}}^j+{\mathrm{\ϵ}}_{\mathrm{cA}}^j---\left(24\right)$

其中，令： $\frac{x_s^j-x_r^0}{{\mathrm{\ρ}}_0^j}=l_j,$ $\frac{y_s^j-y_r^0}{{\mathrm{\ρ}}_0^j}=m_j,$ $\frac{z_s^j-z_r^0}{{\mathrm{\ρ}}_0^j}=n_j$

(4.3)误差方程表示为：

$V_j=-l_j{V}_x-m_j{V}_y-n_j{V}_z-{\mathrm{\λ}}_1{N}_a^j-C({\mathrm{dt}}^j-\mathrm{dT})+R^j---\left(25\right)$

其中： $R^j={\mathrm{\ρ}}_0^j-L_a^j+d_{{\mathrm{trop}}^j}+{\mathrm{drel}}^j$

对于同一历元，对GPS卫星i的观测方程为：

$V_i=-l_i{V}_x-m_i{V}_y-n_i{V}_z-{\mathrm{\λ}}_1{N}_a^i-C({\mathrm{dt}}^i-\mathrm{dT})+R^i---\left(26\right).$

所述的步骤5)中天线相位中心偏差改正，包括如下具体步骤：

星固系的Z轴指向地心，其单位向量e_z为：

$e_z=\frac{-r_{\mathrm{sat}}}{{|r}_{\mathrm{sat}}|}---\left(27\right)$

式中的r_sat为卫星质量中心的坐标；

星固系的Y轴是卫星方向与太阳方向至卫星方向的叉乘，其单位向量e_y为：

$e_y=\frac{r_{\mathrm{sat}}\×(r_{\mathrm{sat}}-r_{\mathrm{sun}})}{|r_{\mathrm{sat}}\×(r_{\mathrm{sat}}-r_{\mathrm{sun}})|}---\left(28\right)$

式中的r_sun为太阳坐标；

星固系的X轴与另外两个轴组成右手系，其单位向量e_x为：

$e_x=\frac{e_y\×e_z}{|e_y\×e_z|}---\left(29\right).$

所述的步骤11)中运用卡尔曼滤波进行参数估计，包括如下具体步骤：

(11.1)假设线性离散系统的状态方程和观测方程为：

            X_k+1＝Φ_k+1，kX_k+Γ_kω_k                (30)

            Z_k+1＝H_k+1X_k+1+v_k+1                    (31)

其中，X_k为n维状态向量；Φ_k+1为n×n维的一步状态转移矩阵；Г_k为n×p维动态噪声驱动阵，ω_k为P维系统动态噪声向量，Z_k+1为m维观测向量，H_k+1为m×n观测矩阵，v_k+1为m维观测噪声向量；

(11.2)首先根据前一次滤波值

(或初值)计算预测值：

$\hat{X}(k/k-1)=\mathrm{\Φ}(k,k-1)\hat{X}(k-1/k-1)---\left(32\right)$

根据前一次得到的滤波误差方差阵P(k-1/k-1)(或初值)及系统噪声的方差阵Q_k计算预测误差方差阵：

P(k/k-1)＝Q(k，k-1)P(k-1/k-1)Q^T(k/k-1)-Г(k-1)Q(k-1)Г^T(k-1)    (33)

(11.3)计算卡尔曼滤波增益，公式如下：

K(k)＝P(k-1/k-1)H^T(k)[H(k)P(k/k-1)H^T(k)+R(k)]^T                  (34)

根据新的观测值z(k)得到：

$v\left(k\right)=z\left(k\right)-H\left(k\right)\hat{X}(k/k-1)---\left(35\right)$

(11.4)计算滤波估计，公式如下：

$\hat{X}(k/k)=\hat{X}(k-1/k-1)+K\left(k\right)v\left(k\right)---\left(36\right)$

计算滤波误差方差阵，公式如下：

                  P(k/k)＝[I-K(k)H(k)]P(k/k-1)                  (37)

(11.5)将滤波估计存入计算机，下一时刻得到新的观测值，重复上述计算过程，滤波以不断“预测-修正”的递推方式进行计算，对预测值不断进行修正。

取2006年1月30日bjfs(北京房山)的L1单差观测值进行试验，所取的卫星号是PRN23和PRN13。bjfs使用的接收机型号为ASHTECH Z-XII3，数据采样率是30s。使用高阶差分的方法对20h 41min 00s到23h 59min 30s共398个星间单差数据进行周跳探测，从图1可以看出，观测值中没有周跳。

通过使用多项式拟合法，得到各个历元的外推值

每次拟合得到的中误差σ(j)(图2)，以及用这些外推值与它们的实际观测值做差得到

(图3)。还可以得到外推值跟实际观测值之差与中误差的比值(图4)

由于星间单差只是消去了接收机的钟差，还剩余卫星钟钟差、电离层、对流层、多路径等因素的影响，以及多项式拟合法还受到采样率的影响。前面选择了k＝5，但当一段观测值序列特别光滑的时候，会出现特别小的情况，此时很难满足公式(3)，而导致似乎存在周跳。为了避免这种情况，根据图3，可以限制σ(j)＞0.2。

根据以上面的数据进行试验。当在23号卫星的L1的第25个历元的载波相位观测值上模拟加上3个周跳时，那么在第24个差值上表现为异常。如图5所示，完全可以通过多项式拟合法探测出来。

然后在两颗卫星的第25个历元的附近用载波相位变化率法进行探测。如图6和图7所示。

可判断出是23号卫星发生了周跳，判断结果达到了预期。

本实施例综合利用多项式拟合法和载波相位变化率法可以准确地确定周跳发生的卫星号和历元，并能修复2周以上的周跳，提高了周跳探测的可靠性和有效性。

实施例2～3：格网模型与半和改正法的解算精度比较

本实施例对2007年1月15日的IGS跟踪站BJFS(bjfs0150.07O)和WUHN(wuhn0150.07O)的数据分别用电离层格网模型和半和改正的方法进行了解算。数据处理中对IGS跟踪站的数据，只利用其观测值中的L1和C/A码数据进行解算。统计结果如图8～11所示。

从图中可以看到：

1)使用格网模型改正，在N、E方向能达到5分米的精度。但是U方向的精度在1米左右；

2)使用半和改正，在N、E、U方向能达到2分米左右的精度；

3)很明显半和改正要比格网模型改正方法所达到的精度要高。其主要原因是半和改正中不需要考虑电离层的影响；而格网模型所改正的电离层延迟并不是很理想，特别是在高程方向的影响还很大。

本发明综合利用了多项式拟合法和载波相位变化率法来探测周跳，结合粗差修改实现单频GPS接收机的精密单点定位。多项式拟合法结合载波相位变化率法可以正确探测与修复周跳。同时，误差改正解算步骤中采用半和改正法，其定位精度在N、E、U方向都能达到2分米左右。上述技术方案有效实现了GPS系统的高精度定位和导航。

Claims (3)

1.一种单频GPS接收机的精密单点定位方法，包括GPS数据的读取和预处理步骤、误差改正解算步骤和结果输出步骤，所述预处理步骤包括进行周跳的探测和修复，其特征在于所述周跳的探测和修复步骤包括：

1)将m个无周跳的载波相位观测值φ_i带入公式(1)，进行多项式拟合；

φ_i＝a₀+a₁(t_i-t₀)+a₂(t_i-t₀)²+...+a_n(t_i-t₀)ⁿi＝1，2，...，m；m＞(n+1)

(1)

用最小二乘法求得式中的多项式系数，并根据拟合后的残差V_i计算出误差：

其中j表示第j次拟合    (2)

2)用求得的多项式系数来外推下一个历元的载波相位观测值

并与实际观测值φ_i进行比较：

其中k为门阀系数    (3)

若公式(3)成立，则认为该观测值没有周跳；

加入上述无周跳的实际观测值后继续上述过程进行多项式拟合；

若公式(4)成立，则认为实际观测值发生了周跳，此时应采用外推的整周计数去取代有周跳的实际观测值中的整周计数，但不足一周的部分F_r(φ)仍保持不变；继续上述过程，直到最后一个观测值为止；

3)利用第i和i-1历元载波相位观测值，获得载波相位观测值的变化率：

φ′_i＝(φ_i-φ_i-1)/Δt    (5)

又由：

φ′_i+1＝φ′_i    (6)

基于第i历元的观测值，并利用上式(6)，可预报第i+1历元的观测值：

若载波相位不出现周跳，它与第i+1历元的观测值的差值应在一定范围之内：

如果公式(8)不成立，则认为有周跳发生并进行修复，精度可在±2周之内；

所述的误差改正解算步骤包括：

4)把L1和C/A码进行组合得到一个新的观测值La，并建立以La为观测值的观测方程，若测站的近似坐标为(x_r ⁰、y_r ⁰、z_r ⁰)，将观测方程在(x_r ⁰、y_r ⁰、z_r ⁰)处用泰勒级数展开得到线性化的观测方程，将误差方程与同一历元中GPS卫星i的观测方程相减，即以GPS卫星i作为参考卫星，得到：

$V_{\mathrm{ij}}=-(l_j-l_i)V_x-(m_j-m_i)V_y-(n_j-n_i)V_z-{\mathrm{\λ}}_1(N_a^j-N_a^i)-C({\mathrm{dt}}^j-{\mathrm{dt}}^i)+(R^j-R^i)---\left(10\right)$

从而消除了电离层一阶项的影响和接收机钟差的影响；

5)卫星质量中心与卫星发射天线的相位中心的不重合产生天线相位中心偏差：该偏差向量在常用的星固坐标系中的三个坐标分量为：星固系的Z轴指向地心，其单位向量为e_z；星固系的Y轴是卫星方向与太阳方向至卫星方向的叉乘，其单位向量为e_y；星固系的x轴与另外两个轴组成右手系，其单位向量为e_x，因此卫星坐标系从质量中心到相位中心的改正值为：

$\mathrm{\Δ}{R}_{\mathrm{sat}}=(e_x,e_y,e_z)\left(\begin{array}{c}\mathrm{\Δ}{X}_{\mathrm{sat}}\\ \mathrm{\Δ}{Y}_{\mathrm{sat}}\\ \mathrm{\Δ}{Z}_{\mathrm{sat}}\end{array}\right)---\left(11\right)$

6)利用伪距做初始解算，获得接收机的近似位置；

7)采用Hopfield模型改正天顶方向对流层延迟：

$\mathrm{\Δ}{P}_d=155.2\×{10}^{-7}\×\frac{P_s}{T_s}\×(h_d-h_s)---\left(12\right)$

$\mathrm{\Δ}{P}_w=155.2\×{10}^{-7}\×\frac{4810}{{T_s}^2}\×e_s\×(h_w-h_s)---\left(13\right)$

h_d＝40136+148.72×(T_s-273.16)    (14)

h_w＝11000    (15)

e_s＝RH·exp(-37.2465+0.21366T-0.0002568908T²)    (16)

其中，T_s，P_s，e_s，h_s分别为测站的温度、气压、湿度、正高，RH为相对湿度，T为干温，为开尔文温度；

8)卫星天线在进出地影时，太阳能面板会慢慢地旋转，导致卫星天线重新定向，改变测站和卫星之间的几何关系，产生相位缠绕，改正方法为：

其中，

是卫星到接收机的单位向量，

是由星固坐标系下的单位向量计算得到的卫星有效的偶极向量，是由地方坐标系下的单位向量

计算得到的接收机有效的偶极向量；

9)利用拉格朗日插值算法对精密星历和精密钟差进行任意时刻的插值：设在n+1个时间为t₁，t₂，...t_n+1插值节点上卫星的坐标分别是x(t₁)，x(t₂)……x(t_n+1)，那么在任意时刻t卫星的坐标用下述公式表示：

$x\left(t\right)=\underset{i=1}{\overset{n+1}{\mathrm{\Σ}}}[x\left(t_i\right)\·\underset{\underset{j\≠i}{j=1}}{\overset{n+1}{\mathrm{\Π}}}\frac{t-t_j}{t_i-t_j}]---\left(21\right)$

利用公式在x、y和z方向上以及根据精密钟差文件分别进行插值处理，即可得到任意时刻卫星的位置和钟差值；

10)万有引力作用下使得地球表面产生固体潮现象，固体潮对测站位置的影响公式如下：

$\mathrm{\Δ}{\stackrel{r}{X}}_i=\underset{i=2}{\overset{3}{\mathrm{\Σ}}}\frac{{\mathrm{GM}}_i}{{\mathrm{GM}}_E}\·\frac{r_E^4}{{\left|X_j\right|}^3}\·\{3\·l_2\·(\frac{X_p}{\left|X_p\right|}\·\frac{X_j}{\left|X_j\right|})\frac{X_j}{\left|X_j\right|}+[3\·(\frac{h_2}{2}-l_2)\·{(\frac{X_p}{\left|X_p\right|}\·\frac{X_j}{\left|X_j\right|})}^2-\frac{h_2}{2}]\·\frac{X_p}{\left|X_p\right|}\}$

$+[-0.025\sin \mathrm{\φ}\cos \mathrm{\φ}\sin (\mathrm{\θ}+\mathrm{\λ})]\·\frac{X_p}{\left|X_p\right|}---\left(22\right)$

其中：r_E为地球的半径；X_j为摄动天体在地心参考框架中的坐标向量；X_p为测站在地心参考框架中的坐标向量；GM_i为摄动天体的引力参数，j＝2时为月球，j＝3时为太阳；GM_E为地球引力参数；h₂为Love数且h₂＝0.6090，l₂为Shida数且l₂＝0.0852；φ，λ为测站纬度和经度；θ为格林尼治平恒星时；

11)运用卡尔曼滤波法，利用受到噪声污染的观测值，对对流层参数、模糊度参数做最优化的估计；

所述的结果输出步骤包括：

12)输出地固系下的(X，Y，Z)或者是大地坐标(B，L，H)，以及在以平均坐标为原点的站心地平坐标下的(N，E，U)计算结果；

所述的步骤4)中利用半合改正来消除电离层影响，包括如下具体步骤：

(4.1)把L1和C/A码进行组合得到一个新的观测值：

其模糊度为

则以La为观测值的观测方程为：

$L_a^j={\mathrm{\ρ}}^j-{\mathrm{\λ}}_1{N}_a^j-C({\mathrm{dt}}^j-\mathrm{dT})+{\mathrm{dtrop}}^j+\mathrm{dmul}{t}_a^j+{\mathrm{drel}}^j+{\mathrm{\ϵ}}_a^j---\left(23\right)$

(4.2)将公式(23)在(x_r ⁰、y_r ⁰、z_r ⁰)处用泰勒级数展开后，得到线性化的观测方程：

$L_a^j={\mathrm{\ρ}}_0^j-\frac{x_s^j-x_r^0}{{\mathrm{\ρ}}_0^j}{V}_x-\frac{y_s^j-y_r^0}{{\mathrm{\ρ}}_0^j}{V}_y-\frac{z_s^j-z_r^0}{{\mathrm{\ρ}}_0^j}{V}_z$

$-{\mathrm{\λ}}_1{N}_a^j-C({\mathrm{dt}}^j-\mathrm{dT})+{\mathrm{dtrop}}^j+{\mathrm{dmult}}_{\mathrm{cA}}^j+{\mathrm{drel}}^j+{\mathrm{\ϵ}}_{\mathrm{cA}}^j---\left(24\right)$

其中，令： $\frac{x_s^j-x_r^0}{{\mathrm{\ρ}}_0^j}=l_j,$ $\frac{y_s^j-y_r^0}{{\mathrm{\ρ}}_0^j}=m_j,$ $\frac{z_s^j-z_r^0}{{\mathrm{\ρ}}_0^j}=n_j$

(4.3)误差方程表示为：

$V_j=-l_j{V}_x-m_j{V}_y-n_j{V}_z-{\mathrm{\λ}}_1{N}_a^j-C({\mathrm{dt}}^j-\mathrm{dT})+R^j---\left(25\right)$

其中： $R^j={\mathrm{\ρ}}_0^j-L_a^j+d_{{\mathrm{trop}}^j}+{\mathrm{drel}}^j$

对于同一历元，对GPS卫星i的观测方程为：

$V_i=-l_i{V}_x-m_i{V}_y-n_i{V}_z-{\mathrm{\λ}}_1{N}_a^i-C({\mathrm{dt}}^i-\mathrm{dT})+R^i---\left(26\right);$

所述的步骤5)中天线相位中心偏差改正，包括如下具体步骤：

星固系的Z轴指向地心，其单位向量e_z为：

$e_z=\frac{-r_{\mathrm{sat}}}{\left|r_{\mathrm{sat}}\right|}---\left(27\right)$

式中的r_sat为卫星质量中心的坐标；

星固系的Y轴是卫星方向与太阳方向至卫星方向的叉乘，其单位向量e_y为：

$e_y=\frac{r_{\mathrm{sat}}\text{\×}(r_{\mathrm{sat}}-r_{\mathrm{sun}})}{|r_{\mathrm{sat}}\×(r_{\mathrm{sat}}-r_{\mathrm{sun}})|}---\left(28\right)$

式中的r_sun为太阳坐标；

星固系的X轴与另外两个轴组成右手系，其单位向量e_x为：

$e_x=\frac{e_y\×e_z}{|e_y\×e_z|}---\left(29\right);$

所述的步骤11)中运用卡尔曼滤波进行参数估计，包括如下具体步骤：

(11.1)假设线性离散系统的状态方程和观测方程为：

X_k+1＝Ф_k+1，kX_k+F_kω_k    (30)

Z_k+1＝H_k+1X_k+1+v_k+1    (31)

其中，X_k为n维状态向量；Ф_k+1为n×n维的一步状态转移矩阵；Γ_k为n×p维动态噪声驱动阵，ω_k为P维系统动态噪声向量，Z_k+1为m维观测向量，H_k+1为m×n观测矩阵，v_k+1为m维观测噪声向量；

(11.2)首先根据前一次滤波值

计算预测值：

$\hat{X}(k/k-1)=\mathrm{\Φ}(k,k-1)\hat{X}(k-1/k-1)---\left(32\right)$

根据前一次得到的滤波误差方差阵P(k-1/k-1)及系统噪声的方差阵Q_k计算预测误差方差阵：

P(k/k-1)＝Q(k，k-1)P(k-1/k-1)Q^T(k/k-1)-Γ(k-1)Q(k-1)Γ^T(k-1)    (33)

(11.3)计算卡尔曼滤波增益，公式如下：

K(k)＝P(k-1/k-1)H^T(k)[H(k)P(k/k-1)H^T(k)+R(k)]^T    (34)

根据新的观测值z(k)得到：

$v\left(k\right)=z\left(k\right)-H\left(k\right)\hat{X}(k/k-1)---\left(35\right)$

(11.4)计算滤波估计，公式如下：

$\hat{X}(k/k)=\hat{X}(k-1/k-1)+K\left(k\right)v\left(k\right)---\left(36\right)$

计算滤波误差方差阵，公式如下：

P(k/k)＝[I-K(k)H(k)]P(k/k-1)    (37)

(11.5)将滤波估计存入计算机，下一时刻得到新的观测值，重复上述计算过程，滤波以不断“预测一修正”的递推方式进行计算，对预测值不断进行修正。

2.如权利要求1所述的单频GPS接收机的精密单点定位方法，其特征在于：所述的步骤1)中，公式(1)中n＝4，m＝14。

3.如权利要求1所述的单频GPS接收机的精密单点定位方法，其特征在于：所述的步骤2)中，公式(3)中k＝5。

Images