CN112147663A - 一种卫星和惯性组合动对动实时精密相对定位方法 - Google Patents

一种卫星和惯性组合动对动实时精密相对定位方法 Download PDF

Info

Publication number
CN112147663A
CN112147663A CN202011327522.4A CN202011327522A CN112147663A CN 112147663 A CN112147663 A CN 112147663A CN 202011327522 A CN202011327522 A CN 202011327522A CN 112147663 A CN112147663 A CN 112147663A
Authority
CN
China
Prior art keywords
satellite
epoch
relative positioning
dynamic
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011327522.4A
Other languages
English (en)
Other versions
CN112147663B (zh
Inventor
李青松
吴杰
董毅
王鼎杰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202011327522.4A priority Critical patent/CN112147663B/zh
Publication of CN112147663A publication Critical patent/CN112147663A/zh
Application granted granted Critical
Publication of CN112147663B publication Critical patent/CN112147663B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/45—Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement
    • G01S19/47—Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement the supplementary measurement being an inertial measurement, e.g. tightly coupled inertial
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/165—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/43—Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry
    • G01S19/44—Carrier phase ambiguity resolution; Floating ambiguity; LAMBDA [Least-squares AMBiguity Decorrelation Adjustment] method

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Automation & Control Theory (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明提出了一种卫星和惯性组合动对动实时精密相对定位方法，基于载波相位实时动态差分相对定位方法获得移动站与动基准之间的相对定位矢量，基于载波相位GNSS/惯导紧组合算法获得测量更新时间间隔△t TC 内的移动站位置增量、动基准位置增量以及移动站在预报时间△t p 内的位置增量，基于滑动窗口的多项式预报方法获得动基准在预报时间，最后合成相对定位结果并输出。该方法不仅可以克服位置增量和GNSS原始观测数据播发延时带来的影响，而且可以通过低数据播发率和采样率提供极高更新率的精密相对定位结果。

Description

一种卫星和惯性组合动对动实时精密相对定位方法
技术领域
本发明属于卫星和惯性组合相对定位导航技术领域，特别涉及一种低数据播发率，低采样率和极高更新率的卫星和惯性组合动对动精密相对定位方法。
背景技术
高更新率的精密相对定位是许多与安全相关的动对动应用中要求的技术，例如车对车协同安全应用、自主空中加油和舰载机着舰等。目前常用的动对动相对定位技术为载波相位RTK（Real-Time Kinematic）差分相对定位技术。然而利用载波相位RTK（Real-TimeKinematic）差分相对定位技术实现高更新率输出主要面临数据播发率大，接收机采样率要求高，播发延时或通信数据包丢失等原因导致卫导观测数据难以同步等问题。
为了克服上述缺点，可以将惯性导航***（INS）测量信息与全球导航卫星***（GNSS）原始观测信息结合来提供高更新率相对定位结果。目前虽已有文献报道卫星/惯性组合相对定位方法（参见[1] Williamson, W. R.; Abdel-Hafez, M. F.; Rhee, I.;Song, E.; Wolfe, J. D.; Chichka, D. F.; Speyer, J. L., An InstrumentationSystem Applied to Formation Flight. IEEE Transactions on Control SystemsTechnology 2007, 15, (1), 75-85. [2] Alam, N.; Kealy, A.; Dempster, A. G., AnINS-Aided Tight Integration Approach for Relative Positioning Enhancement inVANETs. IEEE Transactions on Intelligent Transportation Systems 2013, 14,(4), 1992-1996. [3] Lee, J. Y.; Kim, H. S.; Choi, K. H.; Lim, J.; Chun, S.;Lee, H. K., Adaptive GPS/INS integration for relative navigation. GpsSolutions 2016, 20, (1), 63-75.）。但是这些方法都需要播发原始惯导测量信息，没有考虑数据播发率和采样率带来的通信压力和计算量负担的问题，也没有考虑利用低数据播发率和采样率实现极高更新率厘米级相对定位解的可能性。
发明内容
针对现有技术存在的缺陷，本发明提出了一种卫星和惯性组合动对动实时精密相对定位方法。本发明是一种基于载波相位实时动态（RTK）差分相对定位方法（DGNSS）、载波相位GNSS/惯导紧组合（TC-GNSS/INS）算法和位置增量多项式预报方法的卫星和惯导组合动对动精密相对定位方法。该方法不仅可以克服位置增量和GNSS原始观测数据播发延时带来的影响，而且可以通过低数据播发率和采样率提供极高更新率的精密相对定位结果。
为实现上述技术目的，本发明的技术方案如下：
卫星和惯性组合动对动实时精密相对定位方法，包括以下步骤：
第一步，移动站接收机、动基准接收机同步采样得到原始GNSS观测信息，分别探测与剔除各自采样所得的原始GNSS观测信息中载波相位周跳和伪距粗差，得到预处理后的移动站接收机和动基准接收机的载波相位和伪距观测信息。
第二步，确定载波相位整周模糊度
Figure 754053DEST_PATH_IMAGE001
，利用载波相位实时动态差分相对定位方法相对定位，获得移动站与动基准之间的相对定位矢量dX RB,DGNSS 。
第三步，计算测量更新时间间隔△t TC 内的移动站位置增量
Figure 835142DEST_PATH_IMAGE002
、动基准位置增量
Figure 529559DEST_PATH_IMAGE003
，同时计算移动站在预报时间间隔△t p 内的位置增量
Figure 823138DEST_PATH_IMAGE004
。
第四步，基于滑动窗口的多项式预报方法获得动基准在预报时间间隔△t p 内的位置增量
Figure 535879DEST_PATH_IMAGE005
。
第五步，根据dX RB,DGNSS 、
Figure 725551DEST_PATH_IMAGE006
、
Figure 156533DEST_PATH_IMAGE007
、
Figure 4534DEST_PATH_IMAGE008
和
Figure 837361DEST_PATH_IMAGE009
合成相对定位结果并输出。具体地，合成相对定位结果dX RB,DGNSS/INS 为：
Figure 994673DEST_PATH_IMAGE010
。
进一步地，本发明的第一步中，故障信息包括载波相位周跳和伪距粗差。对于载波相位周跳，采用GF方法与惯导辅助周跳探测相结合来同时探测和排除不同卫星的单频或双频载波相位周跳；对于伪距粗差，采用单站和双站的惯导辅助粗差探测方法探测后剔除。
进一步地，本发明的第二步中，先判断载波相位整周模糊度
Figure 912950DEST_PATH_IMAGE011
是否已知。如果载波相位整周模糊度
Figure 299063DEST_PATH_IMAGE012
已知，则直接利用载波相位实时动态差分相对定位方法相对定位，获得移动站与动基准之间的相对定位矢量dX RB,DGNSS 。如果载波相位整周模糊度
Figure 189659DEST_PATH_IMAGE013
未知，则基于预处理之后的移动站接收机和动基准接收机的载波相位和伪距观测信息，利用本领域公知的最小二乘模糊度降相关平差方法（LAMBDA，参见Teunissen, P. J. G., Theleast-squares ambiguity decorrelation adjustment: a method for fast GPSinteger ambiguity estimation. Journal of Geodesy 1995, 70, (1), 65-82.）求解载波相位整周模糊度
Figure 517872DEST_PATH_IMAGE014
，然后利用载波相位实时动态差分相对定位方法相对定位，获得移动站与动基准之间的相对定位矢量dX RB,DGNSS 。
具体地，移动站与动基准之间的相对定位矢量dX RB,DGNSS 的获取方法如下：
首先，站间星间双差载波相位表示为：
Figure 657867DEST_PATH_IMAGE015
其中DD算子表达式为
Figure 362517DEST_PATH_IMAGE016
，即
Figure 652160DEST_PATH_IMAGE017
，
Figure 151274DEST_PATH_IMAGE018
，
Figure 778565DEST_PATH_IMAGE019
，
Figure 21327DEST_PATH_IMAGE020
；
Figure 433985DEST_PATH_IMAGE021
、
Figure 104001DEST_PATH_IMAGE022
分别表示移动站R观测到的卫星i和j的载波相位观测值；
Figure 421850DEST_PATH_IMAGE023
、
Figure 202724DEST_PATH_IMAGE024
分别表示动基准B观测到的卫星i和j的载波相位观测值；
Figure 984735DEST_PATH_IMAGE025
、
Figure 576385DEST_PATH_IMAGE026
分别表示移动站R观测到的卫星i和j的星地距离；
Figure 178267DEST_PATH_IMAGE027
、
Figure 762832DEST_PATH_IMAGE028
分别表示动基准B观测到的卫星i和j的星地距离；
Figure 133771DEST_PATH_IMAGE029
、
Figure 896322DEST_PATH_IMAGE030
分别表示移动站R观测到的卫星i和j的载波相位模糊度；
Figure 719921DEST_PATH_IMAGE031
、
Figure 311440DEST_PATH_IMAGE032
分别表示动基准B观测到的卫星i和j的载波相位模糊度；
Figure 536884DEST_PATH_IMAGE033
、
Figure 985183DEST_PATH_IMAGE034
分别表示移动站R观测到的卫星i和j的测量误差；
Figure 315320DEST_PATH_IMAGE035
、
Figure 241688DEST_PATH_IMAGE036
分别表示动基准B观测到的卫星i和j的测量误差。
忽略站间视线矢量变化，线性化之后的观测方程表示为：
Figure 321639DEST_PATH_IMAGE037
其中
Figure 144102DEST_PATH_IMAGE038
和
Figure 942294DEST_PATH_IMAGE039
分别表示移动站R观测到的卫星
Figure 157506DEST_PATH_IMAGE040
和
Figure 357543DEST_PATH_IMAGE041
的视线矢量，b RB 表示待求的基线矢量，即移动站与动基准之间的相对定位矢量。
给定至少4颗共视星，移动站与动基准之间的相对位置矢量b RB 由最小二乘方法获得，将获得的移动站与动基准之间的相对位置矢量b RB 记为dX RB,DGNSS 。其中的载波相位模糊度
Figure 944382DEST_PATH_IMAGE042
、
Figure 167553DEST_PATH_IMAGE043
，基于双差载波相位和伪距利用本领域公知的LAMBDA方法（参见Teunissen,P. J. G., The least-squares ambiguity decorrelation adjustment: a method forfast GPS integer ambiguity estimation. Journal of Geodesy 1995, 70, (1), 65-82.）计算获得。
进一步地，本发明中第三步的实现方法如下：
首先确定状态矢量，状态矢量包含9个导航误差状态和6个敏感器偏差状态，表示为：
Figure 452035DEST_PATH_IMAGE044
其中
Figure 506578DEST_PATH_IMAGE045
、
Figure 202002DEST_PATH_IMAGE046
、
Figure 240365DEST_PATH_IMAGE047
、
Figure 62959DEST_PATH_IMAGE048
和
Figure 972009DEST_PATH_IMAGE049
分别表示k历元的位置、速度、姿态误差矢量、加表和陀螺零偏误差。
使用经过第一步故障探测与剔除之后的GNSS原始观测信息构造观测方程，则卫星i在k历元线性化后的观测信息矢量
Figure 103913DEST_PATH_IMAGE050
表示为：
Figure 301676DEST_PATH_IMAGE051
其中r表示参考卫星编号，
Figure 177228DEST_PATH_IMAGE052
表示时间差分算子，
Figure 954167DEST_PATH_IMAGE053
表示线性化后k历元的待探测的卫星i和参考卫星r星间单差伪距，
Figure 256972DEST_PATH_IMAGE054
线性化后k历元的星间单差载波相位的时间差分项。
假设待探测的卫星总数为m，k历元的观测矩阵为：
Figure 4348DEST_PATH_IMAGE055
时间更新和测量更新方程表示为：
Figure 621274DEST_PATH_IMAGE056
其中
Figure 504917DEST_PATH_IMAGE057
表示k-1历元到k历元的状态转移矩阵，W k-1表示k-1历元到k历元的过程噪声矩阵，H k 表示k历元的设计矩阵，V k 表示k历元的测量噪声矩阵，X k-1表示k-1历元的状态矢量。
假设P表示方差矩阵，进行卡尔曼滤波，从k到 k+△t TC 历元，测量更新时间间隔为△t TC 的状态矢量增量
Figure 729356DEST_PATH_IMAGE058
及其相应的方差矩阵
Figure 698449DEST_PATH_IMAGE059
的表达式如下：
Figure 119066DEST_PATH_IMAGE060
Figure 857215DEST_PATH_IMAGE061
其中
Figure 767402DEST_PATH_IMAGE062
Figure 427054DEST_PATH_IMAGE063
Figure 198832DEST_PATH_IMAGE064
Figure 57066DEST_PATH_IMAGE065
Figure 75838DEST_PATH_IMAGE066
Figure 19523DEST_PATH_IMAGE067
表示k历元到k+△t TC 历元的状态转移矩阵，I表示单位矩阵，
Figure 578680DEST_PATH_IMAGE068
表示k+△t TC 历元方差矩阵，Q k-1表示k-1 历元到k历元的过程噪声协方差矩阵，R k 表示k历元的测量噪声协方差矩阵。
利用
Figure 494684DEST_PATH_IMAGE069
和
Figure 231827DEST_PATH_IMAGE070
的表达式分别在移动站和动基准计算得到移动站k到k+△t TC 历元状态矢量增量
Figure 928387DEST_PATH_IMAGE071
，动基准k到k+△t TC 历元状态矢量增量
Figure 228919DEST_PATH_IMAGE072
和移动站k+△t TC 到k+△t TC +△t p 历元状态矢量增量
Figure 796166DEST_PATH_IMAGE073
。其计算表达式为
Figure 219057DEST_PATH_IMAGE074
式中
Figure 613699DEST_PATH_IMAGE075
和
Figure 717921DEST_PATH_IMAGE076
分别表示移动站和动基准k历元到k+△t TC 历元的状态转移矩阵，
Figure 405254DEST_PATH_IMAGE077
表示移动站k+△t TC 到 k+△t TC +△t p 历元的状态转移矩阵，
Figure 999047DEST_PATH_IMAGE078
和
Figure 342303DEST_PATH_IMAGE079
分别表示移动站和动基准k到k+△t TC 历元的过程噪声，
Figure 781375DEST_PATH_IMAGE080
表示移动站k+△t TC 到 k+△t TC +△t p 历元的过程噪声，
Figure 339526DEST_PATH_IMAGE081
、
Figure 41903DEST_PATH_IMAGE082
和
Figure 669194DEST_PATH_IMAGE083
的计算方法与上述
Figure 911956DEST_PATH_IMAGE084
相同，I表示单位矩阵。基于以上计算的状态矢量增量，分别提取
Figure 839461DEST_PATH_IMAGE085
、
Figure 260209DEST_PATH_IMAGE086
和
Figure 312479DEST_PATH_IMAGE087
的前三维构成的向量即为需要获得的矩阵测量更新时间间隔△t TC 内，移动站k到k+△t TC 历元位置增量
Figure 358932DEST_PATH_IMAGE088
，动基准k到k+△t TC 历元位置增量
Figure 140943DEST_PATH_IMAGE089
和移动站在预报时间间隔△t p 内，k+△t TC 到 k+△t TC +△t p 历元位置增量
Figure 185123DEST_PATH_IMAGE090
。
进一步地，本发明的第四步中基于滑动窗口的多项式预报方法中一维方向多项式模型定义为：
Figure 521426DEST_PATH_IMAGE091
其中
Figure 122303DEST_PATH_IMAGE092
表示待建模的变量，T i 表示相对于滑动窗口起始的相对时间， o表示阶数，
Figure 493241DEST_PATH_IMAGE093
、
Figure 708322DEST_PATH_IMAGE094
...
Figure 63080DEST_PATH_IMAGE095
表示多项式系数，
Figure 185757DEST_PATH_IMAGE096
和T i 表达式是：
Figure 880043DEST_PATH_IMAGE097
其中t 1,real 表示滑动窗口第一个位置增量的真实尾部时间，t i,real 表示滑动窗口第i个位置增量的真实尾部时间，
Figure 810566DEST_PATH_IMAGE098
表示t i,real 时刻对应的基准站在测量更新时间间隔△t TC 内的位置增量。
当位置矢量三维的滑动窗口的元素个数大于o+1时, 所有的多项式系数可利用最小二乘方法求解。假设k+△t TC 历元相对于滑动窗口起始的相对时间为t c,real ，则动基准在预报时间间隔△t p 对应的相对时间为t c,real +△t p ，此时动基准位置增量
Figure 387040DEST_PATH_IMAGE099
三维位置均通过下式计算得到：
Figure 313408DEST_PATH_IMAGE100
式中
Figure 862201DEST_PATH_IMAGE101
、
Figure 481401DEST_PATH_IMAGE102
...
Figure 30325DEST_PATH_IMAGE103
、
Figure 963646DEST_PATH_IMAGE104
、
Figure 163684DEST_PATH_IMAGE105
...
Figure 688206DEST_PATH_IMAGE106
和
Figure 239273DEST_PATH_IMAGE107
、
Figure 523755DEST_PATH_IMAGE108
...
Figure 578299DEST_PATH_IMAGE109
分别表示
Figure 476984DEST_PATH_IMAGE110
的x、y和z三维方向的多项式系数，
Figure 249768DEST_PATH_IMAGE111
、
Figure 587209DEST_PATH_IMAGE112
和
Figure 246991DEST_PATH_IMAGE113
分别表示
Figure 582158DEST_PATH_IMAGE114
的x、 y和z三维方向在时间间隔△t p 内的动基准位置增量预报值。
与现有技术相比，本发明具有以下优点：
本发明所提供的方法可以为动对动应用输出与惯性测量信息采样率相等的极高更新率（>100Hz）的厘米级相对位置。
本发明可以通过低数据播发率和采样率提供极高更新率的精密相对定位结果。
本发明可以克服通信延时对实时性能的影响并具有一定容忍数据丢失的能力。
附图说明
图1为本发明的总体框图；
图2为本发明单站故障探测与剔除流程图；
图3为本发明双站故障探测与剔除流程图；
图4为本发明103672.113s时刻相对位置组合原理图；
图5为本发明某次动态试验的水平方向相对定位误差图；
图6为本发明某次动态试验的垂直方向相对定位误差图。
具体实施方式
为了使本发明的技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用于解释本发明，并不用于限定本发明。
实施例1：
参照图1，一种卫星和惯性组合动对动实时精密相对定位方法，包括以下步骤：
第一步，移动站接收机、动基准接收机同步采样得到原始GNSS观测信息，分别探测与剔除各自采样所得的原始GNSS观测信息中载波相位周跳和伪距观测信息的故障信息粗差，得到预处理后的移动站接收机和动基准接收机的载波相位和伪距观测信息。
本发明采用多方法联合故障探测与剔除。GNSS观测数据粗差必须在被用于相对定位之前监测和剔除。故障信息包括载波相位周跳和伪距粗差。载波相位周跳和伪距粗差是最常见的故障。
对于载波相位周跳，经典的GF方法与惯导辅助周跳探测相结合来同时探测和排除不同卫星的单频或双频周跳。GF检验统计量及其方差可表示为
Figure 576659DEST_PATH_IMAGE115
（1）
其中i表示待探测的卫星编号；f 1和f 2表示信号的不同频点，
Figure 717790DEST_PATH_IMAGE116
和
Figure 481347DEST_PATH_IMAGE117
分别表示f 1和 f 2频点的波长；△ t 表示时间差分算子；
Figure 537814DEST_PATH_IMAGE118
和
Figure 222873DEST_PATH_IMAGE119
分别表示卫星
Figure 636537DEST_PATH_IMAGE120
频点f 1和f 2的载波相位观测值，
Figure 520179DEST_PATH_IMAGE121
和
Figure 993886DEST_PATH_IMAGE122
分别表示
Figure 979291DEST_PATH_IMAGE118
和
Figure 196645DEST_PATH_IMAGE119
时间差分后的误差方差。
Figure 138057DEST_PATH_IMAGE118
和
Figure 782665DEST_PATH_IMAGE119
可以是非差，单差或双差形式。非差表示原始测量值不进行组合的运算。单差分为星间单差和站间单差两种形式。星间单差表示将同一时刻不同卫星的载波相位测量值进行差分组合。站间单差表示将同一时刻不同站点接收的载波相位测量值进行差分组合。双差表示在站间单差的基础上，进一步将同一时刻不同卫星的站间单差载波相位测量值进行差分组合。分别将这四种组合形式的载波相位观测值代入式（1）计算检验统计量及其方差，并进一步得到检验门限。将检验统计量与检验门限进行比较：若检验统计量的绝对值小于等于检验门限则检验通过，无故障；若检验统计量的绝对值大于检验门限则检验不通过，有故障。这样就形成了四种形式进行周跳探测的GF方法：非差GF方法，星间单差GF方法，站间单差GF方法和站间星间双差GF方法。GF方法对于探测小周跳具有良好的性能，但是不能探测特殊的双频周跳。此外，该方法也不能区分单频周跳，并且需要双频接收机。
惯导辅助周跳探测方法可以克服GF方法的缺点，该方法的检验统计量及其方差表达式为
Figure 504633DEST_PATH_IMAGE123
（2）
其中r表示参考卫星编号；v表示接收机编号；c表示光速；f表示频点；
Figure 463362DEST_PATH_IMAGE124
表示f频点的波长；
Figure 72329DEST_PATH_IMAGE125
表示卫星i与r星间单差卫星钟差；
Figure 887838DEST_PATH_IMAGE126
表示接收机v中待探测的卫星i与参考卫星r星间单差f频点载波相位观测值；
Figure 34786DEST_PATH_IMAGE127
表示
Figure 593943DEST_PATH_IMAGE128
时间差分后的测量方差；
Figure 572263DEST_PATH_IMAGE129
表示接收机v中卫星i与r星间单差星地距离；
Figure 309406DEST_PATH_IMAGE130
表示
Figure 678071DEST_PATH_IMAGE131
时间差分后的误差方差。
Figure 40919DEST_PATH_IMAGE132
可表示为：
Figure 873746DEST_PATH_IMAGE133
（3）
其中d i 表示卫星i的位置；d r 表示参考卫星r的位置；
Figure 234320DEST_PATH_IMAGE134
表示INS提供的接收机v的位置。该方法利用上述方差计算检验门限，然后将检验统计量与检验门限进行比较：若检验统计量的绝对值小于等于检验门限则检验通过，无故障；若检验统计量的绝对值大于检验门限则检验不通过，有故障。为了表述方便, 由这种检验统计量及其方差构成的周跳探测方法简写为星间单差惯辅周跳探测方法。
当动基准的观测信息和绝对位置被移动站接收时，双站联合的观测信息构成的检验统计量及其方差为
Figure 887018DEST_PATH_IMAGE135
（4）
其中
Figure 801360DEST_PATH_IMAGE136
接收机v和w的卫星i和r之间的站间星间双差星地距离表示为：
Figure 488693DEST_PATH_IMAGE137
其中
Figure 20169DEST_PATH_IMAGE138
表示通过数据链接收到的接收机w的位置，
Figure 425742DEST_PATH_IMAGE139
表示接收机v和w的卫星i和r之间的站间星间双差f频点载波相位观测值，
Figure 864814DEST_PATH_IMAGE140
表示
Figure 609916DEST_PATH_IMAGE141
时间差分后的测量方差，
Figure 859763DEST_PATH_IMAGE142
表示
Figure 752632DEST_PATH_IMAGE143
时间差分后的误差方差，其他定义与前述相同。为了表述方便,基于这种检验统计量及其方差构成的周跳探测方法简写为站间星间双差惯辅周跳探测方法。结合GF方法和惯导辅助方法，各种类型的周跳都可以被探测和剔除。
对于伪距粗差，考虑到纯卫导方法的性能较差，本实施例只采用单站和双站的惯导辅助粗差探测方法。单站检验统计量及其方差为
Figure 198657DEST_PATH_IMAGE144
（5）
其中i表示待探测的卫星编号，r表示参考卫星编号，v表示接收机编号；c表示光速，f表示信号频点，
Figure 595004DEST_PATH_IMAGE145
表示接收机
Figure 530599DEST_PATH_IMAGE146
中卫星i与 r星间单差f频点伪距观测值，
Figure 395917DEST_PATH_IMAGE147
表示卫星i与 r星间单差卫星钟差，
Figure 645633DEST_PATH_IMAGE148
表示表示接收机v中卫星i与 r星间单差星地距离，与前述定义相同，
Figure 162065DEST_PATH_IMAGE149
表示接收机v中卫星i与 r星间单差f频点的电离层延迟，利用双差消电离层组合消除,
Figure 2982DEST_PATH_IMAGE150
表示接收机v中卫星i与r星间单差f频点的对流层延迟，使用Saastamoinen模型补偿，表示
Figure 392692DEST_PATH_IMAGE152
时间差分后的误差方差，与前述定义相同，
Figure 514363DEST_PATH_IMAGE153
表示
Figure 791761DEST_PATH_IMAGE154
时间差分后的测量方差。为了表述方便, 基于这种检验统计量及其方差构成的周跳探测方法简写为星间单差惯辅粗差探测方法。
在短基线条件下，电离层和对流层延迟可假设通过双差组合消除，从而双站检验统计量及其方差为
Figure 880940DEST_PATH_IMAGE155
（6）
式中
Figure 3616DEST_PATH_IMAGE156
表示表示接收机v和w的卫星i和 r之间的站间星间双差f频点伪距观测值，
Figure 432324DEST_PATH_IMAGE157
表示接收机v和w的卫星i和 r之间的站间星间双差星地距离，与前述定义相同，
Figure 646003DEST_PATH_IMAGE158
表示
Figure 222478DEST_PATH_IMAGE159
误差方差，
Figure 352108DEST_PATH_IMAGE160
表示
Figure 697639DEST_PATH_IMAGE161
误差方差，与前述定义相同。为了表述方便, 基于这种检验统计量及其方差构成的周跳探测方法简写为站间星间双差惯辅粗差探测方法。
以上所有检验统计量在原始假设中都服从方差已知的0均值整体分布。在备择假设中，均值等于粗差并且方差保持不变。因此，给定要求的误警率P FA ，可求得检验门限为
Figure 316839DEST_PATH_IMAGE162
（7）
其中σ T 表示检验统计量标准差Φ-1(x) 表示 Φ(x)逆函数, 其定义为
Figure 865763DEST_PATH_IMAGE163
（8）
可以将上述所有方法结合起来探测单站和双站组合观测信息的载波相位周跳和伪距粗差。
为了充分发挥各方法的优点，本实施例分别为单站和双站设计了一种联合的故障探测规则，如图2 和3 所示。
图2中单站联合故障探测步骤为：
①使用非差GF方法初步探测载波相位观测信息的周跳。
②在非差GF方法检测后无周跳的卫星中选择仰角最大的卫星作为参考星。
③ 基于惯导***预报的位置信息利用星间单差惯辅周跳探测方法探测载波相位观测信息的周跳。
④ 利用星间单差GF方法探测载波相位信息周跳。
⑤ 综合星间单差惯辅周跳探测方法和星间单差GF方法周跳探测结果。判断两种方法中检验是否都通过：若是，则无周跳；若否，则有周跳，将该卫星观测信息标注为无效。
⑥ 基于惯导***预报的位置信息利用星间单差惯辅粗差探测方法探测伪距观测信息的粗差。判断检验是否通过：若是，则无粗差；若否，则有粗差，将该卫星观测信息标注为无效。
图3中双站联合故障探测步骤为：
①使用站间单差GF方法初步探测载波相位观测信息的周跳。
②在站间单差GF方法检测后无周跳的卫星中选择仰角最大的卫星作为参考星。
③ 基于惯导***预报的位置信息利用站间星间双差惯辅周跳探测方法探测载波相位观测信息的周跳。
④ 利用站间星间双差GF方法探测载波相位信息周跳。
⑤ 综合站间星间双差惯辅周跳探测方法和站间星间双差GF方法周跳探测结果。判断两种方法中检验是否都通过：若是，则无周跳；若否，则有周跳，将该卫星观测信息标注为无效。
⑥ 基于惯导***预报的位置信息利用站间星间双差惯辅粗差探测方法探测伪距观测信息的粗差。判断检验是否通过：若是，则无粗差；若否，则有粗差，将该卫星观测信息标注为无效。
第二步，判断载波相位整周模糊度
Figure 799084DEST_PATH_IMAGE164
是否已知。若否，基于预处理之后的移动站接收机和动基准接收机的载波相位和伪距观测信息，利用本领域公知的最小二乘模糊度降相关平差方法（LAMBDA）求解载波相位整周模糊度
Figure 999121DEST_PATH_IMAGE164
，然后利用载波相位实时动态差分相对定位方法相对定位，获得移动站与动基准之间的相对定位矢量dX RB,DGNSS ；若是，直接利用载波相位实时动态差分相对定位方法相对定位，获得移动站与动基准之间的相对定位矢量dX RB,DGNSS 。
由于电离层和对流层误差在移动站与动基准之间短基线（< 10 km）高度相关，利用站间星间双差载波相位可以消除接收机钟差，卫星钟差，电离层延迟和对流层延迟等其他公共误差。
站间星间双差载波相位表示为：
Figure 789223DEST_PATH_IMAGE165
（9）
其中DD算子表达式为
Figure 74711DEST_PATH_IMAGE166
，即
Figure 546143DEST_PATH_IMAGE167
，
Figure 351419DEST_PATH_IMAGE168
，
Figure 578001DEST_PATH_IMAGE169
，
Figure 350785DEST_PATH_IMAGE170
；φ表示载波相位观测值，ρ表示星地距离，λ表示信号波长，N表示载波相位模糊度，ε表示测量误差，下标R和B分别表示移动站和动基准，上标i和j表示待探测的卫星编号。
忽略站间视线矢量变化，线性化之后的观测方程可表示为
Figure 360330DEST_PATH_IMAGE171
（10）
其中
Figure 534959DEST_PATH_IMAGE172
和
Figure 417596DEST_PATH_IMAGE173
分别表示卫星i和j 的视线矢量，b RB 表示待求的基线矢量，即移动站与动基准之间的相对定位矢量。载波相位模糊度可以基于双差载波相位和伪距利用本领域公知的LAMBDA方法（参见Teunissen, P. J. G., The least-squares ambiguitydecorrelation adjustment: a method for fast GPS integer ambiguity estimation.Journal of Geodesy 1995, 70, (1), 65-82.）计算获得。
给定至少4颗共视星，移动站与动基准之间的相对位置矢量b RB 可由最小二乘方法（LSE）获得，将获得的移动站与动基准之间的相对位置矢量b RB 记为dX RB,DGNSS 。
第三步，利用TC-GNSS/INS算法计算测量更新时间间隔
Figure 677676DEST_PATH_IMAGE174
内的移动站位置增量
Figure 490911DEST_PATH_IMAGE175
和动基准位置增量
Figure 520047DEST_PATH_IMAGE176
。同时利用TC-GNSS/INS算法计算移动站在预报时间间隔
Figure 88431DEST_PATH_IMAGE177
内的位置增量
Figure 318031DEST_PATH_IMAGE178
。
首先确定状态矢量。
为了提高位置增量的精度，在测量方程中采用载波相位和伪距观测值。状态矢量包含9个导航误差状态和6个敏感器偏差状态，可表示为
Figure 731695DEST_PATH_IMAGE179
（11）
其中
Figure 818599DEST_PATH_IMAGE180
、
Figure 292306DEST_PATH_IMAGE181
、
Figure 261399DEST_PATH_IMAGE182
、
Figure 963907DEST_PATH_IMAGE183
和ε k 分别表示k历元的位置、速度、姿态误差矢量、加表和陀螺零偏误差。
使用的经过第一步故障探测与剔除之后的GNSS原始观测信息构造观测方程，则卫星i在k历元线性化后的观测信息矢量
Figure 233214DEST_PATH_IMAGE184
可表示为：
Figure 877822DEST_PATH_IMAGE185
（12）
其中r表示参考卫星编号，△ t 表示时间差分算子，
Figure 537474DEST_PATH_IMAGE186
表示线性化后k历元的待探测的卫星i和参考卫星r星间单差伪距，
Figure 292940DEST_PATH_IMAGE187
线性化后k历元的星间单差载波相位的时间差分项。
对于载波相位，历元间差分用于去除模糊度参数。对于伪距，电离层延迟利用双频消电离层组合消除，对流层延迟利用Saastamoinen 模型补偿。
假设待探测的卫星总数为m，k历元的观测矩阵为：
Figure 901907DEST_PATH_IMAGE188
（13）
时间更新和测量更新方程表示为：
Figure 717416DEST_PATH_IMAGE189
（14）
其中
Figure 864364DEST_PATH_IMAGE190
表示k-1历元到k历元的状态转移矩阵，W k-1表示k-1历元到k历元的过程噪声矩阵，H k 表示k历元的设计矩阵，V k 表示k历元的测量噪声矩阵，X k-1表示k-1历元的状态矢量。
假设P表示方差矩阵，卡尔曼滤波算法可以进行。从k到 k+△t TC 历元，测量更新时间间隔为△t TC 的状态矢量增量
Figure 423521DEST_PATH_IMAGE191
及其相应的方差矩阵
Figure 401842DEST_PATH_IMAGE192
的表达式如下：
Figure 873405DEST_PATH_IMAGE193
（15）
Figure 507649DEST_PATH_IMAGE194
（16）
其中
Figure 870497DEST_PATH_IMAGE195
（17）
Figure 703324DEST_PATH_IMAGE196
（18）
Figure 595057DEST_PATH_IMAGE197
（19）
Figure 716596DEST_PATH_IMAGE198
（20）
Figure 371218DEST_PATH_IMAGE199
（22）
式中
Figure 58552DEST_PATH_IMAGE200
表示k历元到k+△t TC 历元的状态转移矩阵，I表示单位矩阵，
Figure 386765DEST_PATH_IMAGE201
表示k+△t TC 历元方差矩阵，Q k-1 表示k-1 历元到k历元的过程噪声协方差矩阵，R k 表示k历元的测量噪声协方差矩阵。
利用
Figure 730022DEST_PATH_IMAGE202
和
Figure 700252DEST_PATH_IMAGE203
的表达式分别在移动站和动基准计算得到移动站k到k+△t TC 历元状态矢量增量
Figure 992824DEST_PATH_IMAGE204
，动基准k到k+△t TC 历元状态矢量增量
Figure 695201DEST_PATH_IMAGE205
和移动站k+△t TC 到k+△t TC +△t p 历元状态矢量增量
Figure 322491DEST_PATH_IMAGE206
。其计算表达式为：
Figure 565253DEST_PATH_IMAGE207
（23）
式中
Figure 227179DEST_PATH_IMAGE208
和
Figure 100457DEST_PATH_IMAGE209
分别表示移动站和动基准k历元到k+△t TC 历元的状态转移矩阵，
Figure 965776DEST_PATH_IMAGE210
表示移动站k+△t TC 到k+△t TC +△t p 历元的状态转移矩阵，W R,K 和 W B,K 分别表示移动站和动基准k到k+△t TC 历元的过程噪声，
Figure 12229DEST_PATH_IMAGE211
表示移动站k+△t TC 到k+△t TC +△t p 历元的过程噪声，
Figure 731924DEST_PATH_IMAGE212
、
Figure 572841DEST_PATH_IMAGE213
和
Figure 440303DEST_PATH_IMAGE214
的计算方法与上述
Figure 775600DEST_PATH_IMAGE215
相同，I表示单位矩阵。
基于以上计算的状态矢量增量
Figure 349801DEST_PATH_IMAGE216
、
Figure 361619DEST_PATH_IMAGE217
和
Figure 450798DEST_PATH_IMAGE218
，分别提取
Figure 839054DEST_PATH_IMAGE219
、
Figure 267761DEST_PATH_IMAGE220
和
Figure 198284DEST_PATH_IMAGE221
的前三维构成的向量即为需要获得的
矩阵测量更新时间间隔△t TC 内，移动站k到k+△t TC 历元位置增量
Figure 40338DEST_PATH_IMAGE222
，动基准k到k+△t TC 历元位置增量
Figure 966706DEST_PATH_IMAGE223
和移动站在预报时间间隔
Figure 249919DEST_PATH_IMAGE224
内，k+△t TC 到k+△t TC +△t p 历元位置增量
Figure 134699DEST_PATH_IMAGE225
。
第四步，基于滑动窗口的多项式预报方法获得动基准在预报时间△t p 内的位置增量
Figure 683623DEST_PATH_IMAGE226
。
为了获得动基准播发间隙同步相对位置，采用了基于滑动窗口的多项式预报方法。
多项式预报方法中一维方向多项式模型定义为：
Figure 616944DEST_PATH_IMAGE227
其中
Figure 816981DEST_PATH_IMAGE228
表示待建模的变量；T i 表示相对于滑动窗口起始的相对时间；o表示阶数；
Figure 607082DEST_PATH_IMAGE229
、
Figure 95833DEST_PATH_IMAGE230
...
Figure 629582DEST_PATH_IMAGE231
表示表示多项式系数。
Figure 434858DEST_PATH_IMAGE232
和
Figure 599123DEST_PATH_IMAGE233
表达式是：
Figure 371907DEST_PATH_IMAGE234
（24）
其中t 1,real 表示滑动窗口第一个位置增量的真实尾部时间，t i,real 表示滑动窗口第i个位置增量的真实尾部时间，
Figure 443769DEST_PATH_IMAGE235
表示t i,real 时刻对应的基准站在测量更新时间间隔△t TC 内的位置增量。
当位置矢量三维的滑动窗口的元素个数大于o+1时, 所有的多项式系数可利用最小二乘方法求解。假设k+△t TC 历元相对于滑动窗口起始的相对时间为t c,real ，则动基准在预报时间间隔△t p 对应的相对时间为t c,real +△t p ，此时动基准位置增量
Figure 821660DEST_PATH_IMAGE236
三维位置均通过下式计算得到：
Figure 953564DEST_PATH_IMAGE237
（25）
式中
Figure 964377DEST_PATH_IMAGE238
、
Figure 777612DEST_PATH_IMAGE239
、
Figure 806748DEST_PATH_IMAGE240
、
Figure 109553DEST_PATH_IMAGE241
...
Figure 60192DEST_PATH_IMAGE242
和
Figure 739435DEST_PATH_IMAGE243
、
Figure 365020DEST_PATH_IMAGE244
...
Figure 838727DEST_PATH_IMAGE245
分别表示
Figure 11082DEST_PATH_IMAGE246
的x、 y和z三维方向的多项式系数，
Figure 228437DEST_PATH_IMAGE247
、
Figure 232165DEST_PATH_IMAGE248
和
Figure 814456DEST_PATH_IMAGE249
分别表示
Figure 287157DEST_PATH_IMAGE250
的x、 y和z三维方向时间间隔△t p 内的动基准位置增量预报值。
第五步，合成相对定位结果dX RB,DGNSS/INS ，并输出，其中：
Figure 42624DEST_PATH_IMAGE251
（26）
实施例2：
为了测试实施例1中的方法，本实施例开展了车对车的实测试验, 移动站和动基准都是车并且在广场上运动时有许多转弯。在试验地点附近已知位置架设了一个静态基准站来分别计算移动站与动基准的事后处理的相对位置。相应的结果用来为位置增量和相对位置提供参考结果。在试验中，由Sensonor STIM300 MEMS 和ComNavOEM-K508板组成的GNSS/MEMS原型***固连在移动站二号动基准。GNSS接收机可以为实时和事后导航提供5个频点的观测信息 (B1/B2/B3/L1/L2) 。后续分析中实际只用到四个频点的观测信息 (B1/B3 /L1/L2)。GNSS接收机最大采样率为10 Hz，MEMS采样率为125Hz. Digi International Inc.生产的Xtend-PKG 900 MHz RF电台被用于播发和接收数据包。GNSS原始观测信息以 RTCM3.2 中的 MSM 5格式播发差分数据信息。为了分析位置增量和组合的同步相对位置的精度，可靠事后处理的移动站与静基准，动基准与静基准 DGNSS相对位置被用来作为参考。试验过程中移动站和动基准的GNSS接收机采样率均设置为1Hz，动基准GNSS原始观测数据播发率为1Hz，动基准位置增量播发率为10Hz。
以移动站MEMS观测时为103672.113s相对定位计算为例。此时处于移动站GNSS观测数据采样间隙，同时又处于动基准GNSS原始观测数据和位置增量播发间隙，因此需要将最近时刻的RTK DGNSS相对位置结果，紧组合外推得到的位置增量和多项式预报得到的位置增量组合获得当前的动对动相对定位结果，具体组合原理如图4所示。
（1）首先计算103672.113s的dX RB,DGNSS 、
Figure 104120DEST_PATH_IMAGE252
和
Figure 185209DEST_PATH_IMAGE253
时用到了移动站接收机采样的103672.0时刻GNSS原始观测信息，动基准接收机采样的103672.0时刻GNSS原始观测信息。在将GNSS观测信息用于定位之前，采用单站和双站相结合的多方法联合故障探测与剔除方法探测载波相位周跳和伪距粗差。经比较，各检验统计量均为超过检验门限，无故障卫星。
（2）此时载波相位整周模糊度已知，基于接收到的103672.0s时刻动基准和移动站的GNSS原始观测信息，利用RTK DGNSS 相对定位方法计算103672.0s时刻同步相对位置dX RB,DGNSS 。此时北斗共视星有7个，GPS共视星有9个，除去两颗参考星，此时双频载波双差观测值有28个，计算可得
Figure 128894DEST_PATH_IMAGE254
。
（3） 利用TC-GNSS/INS算法计算动基准和移动站位置增量。
动基准利用TC-GNSS/INS算法计算
Figure 625735DEST_PATH_IMAGE255
，测量更新频率为1Hz，外推位置增量以10Hz保存并播发，此时移动站需要使用接收到的103672.1时刻的位置增量
Figure 89208DEST_PATH_IMAGE256
，此时
Figure 341198DEST_PATH_IMAGE257
。同时由于移动站测量更新频率为1 Hz，而103672.113s离整秒时刻很近，因此可利用TC-GNSS/INS算法直接从103672.0s时刻外推到当前，此时
Figure 772179DEST_PATH_IMAGE258
即此时直接得到
Figure 72711DEST_PATH_IMAGE259
。由计算过程可知
Figure 639958DEST_PATH_IMAGE260
。
（4） 动基准位置增量多项式预报。
由于此时处于动基准位置增量播发间隙，因此需要采用多项式预报方法计算
Figure 813582DEST_PATH_IMAGE261
。多项式算法中阶数取为二阶，滑动窗口历元个数取为10。基于滑动窗口数据，利用最小二乘法可得各向多项式系数为：
Figure 935121DEST_PATH_IMAGE262
当前时刻相对于滑动窗口第一个位置增量的相对时间为1.013s，依据多项式模型可得三个时刻多项式建模参数值为：
Figure 570502DEST_PATH_IMAGE263
又预报时间
Figure 523415DEST_PATH_IMAGE264
，从而可得三个方向预报的位置增量为：
Figure 54890DEST_PATH_IMAGE265
（5）极高更新率相对定位结果合成。
基于前述计算的各分量结果，可以合成当前相对位置为：
Figure 460464DEST_PATH_IMAGE266
上述时刻算例直观地展示了本发明所提供的方法实现极高更新率相对定位的过程。
图5和图6 分别给出了水平方向和垂直方向相对定位误差。图5给出了水平方向相对定位误差散点图，可以发现所有的点都在0.1m散布以内，大部分点散布在0.05m以内。图6给出了垂直方向误差随试验时间变化图，可知垂直方向误差绝大部分小于0.1m。从图中可以看出，在低数据播发率和采样率条件下，相对定位误差仍可以保持在厘米级。
表1 所选数据播发率和采样率条件下相对定位精度
Figure 912917DEST_PATH_IMAGE267
表1显示了相对定位精度统计结果。可以看出无论是各向定位精度还是三维定位精度都非常高，进一步证明了本发明方法的高精度性能。
表2 本发明方法与10Hz纯RTK DGNSS相对定位方法对比结果
Figure 658019DEST_PATH_IMAGE268
表2 显示了本发明方法与10Hz纯RTK DGNSS相对定位方法对比结果。其中，常用的10Hz纯RTK DGNSS相对定位方法简记为方法1，本发明方法简记为方法2。由表可知，本发明方法的数据播发率和接收机采样率都有大幅减少，分别减少了76.24%和90%。而输出更新率有了大幅增加，增加了12.5倍。
表3 不同位置增量播发最小延时条件下统计精度
Figure 422713DEST_PATH_IMAGE269
表4 不同GNSS原始观测信息播发最小延时条件下统计精度
Figure 50003DEST_PATH_IMAGE270
为了验证本发明方法的实时性能，根据记录的数据事后分析中仿真了位置增量和GNSS原始观测信息时延来分析。位置增量和GNSS原始观测信息都设置了最小时延。
表3和4给出了位置增量和GNSS原始观测信息不同最小时延条件下的统计精度。显然相对位置精度会随着仿真的最小时延的增加而变差。对于位置增量，当最小时延小于或等于0.4s时，各向精度可以保持在厘米级。这意味着模式3中所提方法可以忍受连续4个位置增量数据包丢失。对于GNSS原始观测信息，不向垂向水平方向衰减很慢。当最小时延小于或等于2.5s时，三个方向精度仍然可以保持在厘米级，这意味着所提方法可以忍受连续2个GNSS观测数据数据包丢失。通常而言通信延迟一般在毫秒量级，有时可能达到100 ms。因此，本发明方法可以有效克服时延带来的影响并且具有忍受数据包丢失的能力。
综上所述，虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明，任何本领域普通技术人员，在不脱离本发明的精神和范围内，当可作各种更动与润饰，因此本发明的保护范围当视权利要求书界定的范围为准。

Claims (10)

1.卫星和惯性组合动对动实时精密相对定位方法，其特征在于，包括以下步骤：
第一步，移动站接收机、动基准接收机同步采样得到原始GNSS观测信息，分别探测与剔除各自采样所得的原始GNSS观测信息中载波相位周跳和伪距粗差，得到预处理后的移动站接收机和动基准接收机的载波相位和伪距观测信息；
第二步，确定载波相位整周模糊度
Figure 816296DEST_PATH_IMAGE001
，利用载波相位实时动态差分相对定位方法相对定位，获得移动站与动基准之间的相对定位矢量dX RB,DGNSS ；
第三步，计算测量更新时间间隔△t TC 内的移动站位置增量
Figure 264595DEST_PATH_IMAGE002
、动基准位置增量
Figure 575490DEST_PATH_IMAGE003
，同时计算移动站在预报时间间隔△t p 内的位置增量
Figure 236279DEST_PATH_IMAGE004
；
第四步，基于滑动窗口的多项式预报方法获得动基准在预报时间间隔△t p 内的位置增量
Figure 69892DEST_PATH_IMAGE005
；
第五步，根据dX RB,DGNSS 、、
Figure 221705DEST_PATH_IMAGE007
、
Figure 686184DEST_PATH_IMAGE008
和
Figure 620642DEST_PATH_IMAGE009
合成相对定位结果并输出。
2.根据权利要求1所述的卫星和惯性组合动对动实时精密相对定位方法，其特征在于：第一步中，故障信息包括载波相位周跳和伪距粗差。
3.根据权利要求2所述的卫星和惯性组合动对动实时精密相对定位方法，其特征在于：第一步中，对于载波相位周跳，采用GF方法与惯导辅助周跳探测相结合来同时探测和排除不同卫星的单频或双频载波相位周跳；对于伪距粗差，采用单站和双站的惯导辅助粗差探测方法探测后剔除。
4.根据权利要求1、2或3所述的卫星和惯性组合动对动实时精密相对定位方法，其特征在于：第二步中，先判断载波相位整周模糊度
Figure 879585DEST_PATH_IMAGE010
是否已知，如果载波相位整周模糊度
Figure 915806DEST_PATH_IMAGE010
已知，则直接利用载波相位实时动态差分相对定位方法相对定位，获得移动站与动基准之间的相对定位矢量dX RB,DGNSS ；如果载波相位整周模糊度
Figure 918397DEST_PATH_IMAGE010
未知，则基于预处理之后的移动站接收机和动基准接收机的载波相位和伪距观测信息，利用最小二乘模糊度降相关平差方法求解载波相位整周模糊度，然后利用载波相位实时动态差分相对定位方法相对定位，获得移动站与动基准之间的相对定位矢量dX RB,DGNSS 。
5.根据权利要求4所述的卫星和惯性组合动对动实时精密相对定位方法，其特征在于：移动站与动基准之间的相对定位矢量dX RB,DGNSS 的获取方法，如下：
站间星间双差载波相位表示为：
Figure 972940DEST_PATH_IMAGE011
其中DD算子表达式为
Figure 668364DEST_PATH_IMAGE012
，即
Figure 441148DEST_PATH_IMAGE013
，，
Figure 907212DEST_PATH_IMAGE015
，
Figure 507958DEST_PATH_IMAGE016
；
Figure 768038DEST_PATH_IMAGE017
、
Figure 378011DEST_PATH_IMAGE018
分别表示移动站R观测到的卫星i和j的载波相位观测值；
Figure 892300DEST_PATH_IMAGE019
、
Figure 929526DEST_PATH_IMAGE020
分别表示动基准B观测到的卫星i和j的载波相位观测值；
Figure 411323DEST_PATH_IMAGE021
、
Figure 559408DEST_PATH_IMAGE022
分别表示移动站R观测到的卫星i和j的星地距离；
Figure 443050DEST_PATH_IMAGE023
、
Figure 651178DEST_PATH_IMAGE024
分别表示动基准B观测到的卫星i和j的星地距离；
Figure 368073DEST_PATH_IMAGE025
、
Figure 54270DEST_PATH_IMAGE026
分别表示移动站R观测到的卫星i和j的载波相位模糊度；
Figure 792418DEST_PATH_IMAGE027
、
Figure 171447DEST_PATH_IMAGE028
分别表示动基准B观测到的卫星i和j的载波相位模糊度；
Figure 627836DEST_PATH_IMAGE029
、
Figure 134035DEST_PATH_IMAGE030
分别表示移动站R观测到的卫星i和j的测量误差；
Figure 461111DEST_PATH_IMAGE031
、
Figure 11041DEST_PATH_IMAGE032
分别表示动基准B观测到的卫星i和j的测量误差，
Figure 954727DEST_PATH_IMAGE033
表示信号波长；
忽略站间视线矢量变化，线性化之后的观测方程表示为：
Figure 248305DEST_PATH_IMAGE034
其中
Figure 961046DEST_PATH_IMAGE035
和
Figure 167030DEST_PATH_IMAGE036
分别表示移动站R观测到的卫星
Figure 598012DEST_PATH_IMAGE037
和
Figure 429701DEST_PATH_IMAGE038
的视线矢量，b RB 表示待求的基线矢量，即移动站与动基准之间的相对定位矢量；
给定至少4颗共视星，移动站与动基准之间的相对位置矢量b RB 由最小二乘方法获得，将获得的移动站与动基准之间的相对位置矢量b RB 记为dX RB,DGNSS 。
6.根据权利要求5所述的卫星和惯性组合动对动实时精密相对定位方法，其特征在于：载波相位模糊度基于双差载波相位和伪距利用最小二乘模糊度降相关平差方法方法计算获得。
7.根据权利要求5所述的卫星和惯性组合动对动实时精密相对定位方法，其特征在于：第三步的实现方法如下：
确定状态矢量，状态矢量包含9个导航误差状态和6个敏感器偏差状态，表示为：
Figure 262528DEST_PATH_IMAGE039
其中
Figure 888682DEST_PATH_IMAGE040
、
Figure 557691DEST_PATH_IMAGE041
、
Figure 193072DEST_PATH_IMAGE042
、
Figure 349247DEST_PATH_IMAGE043
和ε k 分别表示k历元的位置、速度、姿态误差矢量、加表和陀螺零偏误差；
使用经过第一步故障探测与剔除之后的GNSS原始观测信息构造观测方程，则卫星i在k历元线性化后的观测信息矢量表示为：
Figure 677460DEST_PATH_IMAGE044
其中r表示参考卫星编号，表示时间差分算子，
Figure 551875DEST_PATH_IMAGE045
表示线性化后k历元的待探测的卫星i和参考卫星r星间单差伪距，
Figure 256526DEST_PATH_IMAGE046
线性化后k历元的星间单差载波相位的时间差分项；
假设待探测的卫星总数为m，k历元的观测矩阵为：
Figure 376247DEST_PATH_IMAGE047
时间更新和测量更新方程表示为：
Figure 609782DEST_PATH_IMAGE048
其中
Figure 237072DEST_PATH_IMAGE049
表示k-1历元到k历元的状态转移矩阵，W k-1表示k-1历元到k历元的过程噪声矩阵，H k 表示k历元的设计矩阵，V k 表示k历元的测量噪声矩阵，X k-1表示k-1历元的状态矢量；
假设P表示方差矩阵，进行卡尔曼滤波，从k到 k+△t TC 历元，测量更新时间间隔为△t TC 的状态矢量增量
Figure 214256DEST_PATH_IMAGE050
及其相应的方差矩阵
Figure 610602DEST_PATH_IMAGE051
的表达式如下：
Figure 765771DEST_PATH_IMAGE052
Figure 880357DEST_PATH_IMAGE053
其中：
Figure 395652DEST_PATH_IMAGE054
Figure 646505DEST_PATH_IMAGE055
Figure 487422DEST_PATH_IMAGE056
Figure 823726DEST_PATH_IMAGE057
Figure 159023DEST_PATH_IMAGE058
Figure 529962DEST_PATH_IMAGE059
表示k历元到k+△t TC 历元的状态转移矩阵，I表示单位矩阵，
Figure 276201DEST_PATH_IMAGE060
表示k+△t TC 历元方差矩阵，Q k-1表示k-1 历元到k历元的过程噪声协方差矩阵，R k 表示k历元的测量噪声协方差矩阵；
利用
Figure 834221DEST_PATH_IMAGE061
和
Figure 222477DEST_PATH_IMAGE062
的表达式分别在移动站和动基准计算得到移动站k到k+△t TC 历元状态矢量增量
Figure 933075DEST_PATH_IMAGE063
，动基准k到k+△t TC 历元状态矢量增量
Figure 850216DEST_PATH_IMAGE064
和移动站k+△t TC 到k+△t TC +△t p 历元状态位置矢量增量
Figure 426690DEST_PATH_IMAGE065
，分别提取
Figure 87479DEST_PATH_IMAGE066
、
Figure 167430DEST_PATH_IMAGE067
和
Figure 521051DEST_PATH_IMAGE068
的前三维构成的向量即为需要获得的测量更新时间间隔△t TC 内，移动站k到k+△t TC 历元位置增量
Figure 67046DEST_PATH_IMAGE069
、动基准k到k+△t TC 历元位置增量
Figure 531525DEST_PATH_IMAGE070
和移动站在预报时间间隔△t p 内，k+△t TC 到k+△t TC +△t p 历元位置增量
Figure 465983DEST_PATH_IMAGE071
。
8.根据权利要求7所述的卫星和惯性组合动对动实时精密相对定位方法，其特征在于：第四步中基于滑动窗口的多项式预报方法中多项式模型定义为：
Figure 990505DEST_PATH_IMAGE072
其中
Figure 10414DEST_PATH_IMAGE073
表示待建模的变量，T i 表示相对于滑动窗口起始的相对时间，o表示阶数，
Figure 29317DEST_PATH_IMAGE074
表示多项式系数，
Figure 818281DEST_PATH_IMAGE075
和T i 表达式是：
Figure 779284DEST_PATH_IMAGE076
其中t 1,real 表示滑动窗口第一个位置增量的真实尾部时间，t i,real 表示滑动窗口第i个位置增量的真实尾部时间，
Figure 286489DEST_PATH_IMAGE077
表示t i,real 时刻对应的基准站在测量更新时间间隔△t TC 内的位置增量。
9.根据权利要求8所述的卫星和惯性组合动对动实时精密相对定位方法，其特征在于：第四步中，当位置矢量三维的滑动窗口的元素个数大于o+1时，所有的多项式系数利用最小二乘方法求解，假设k+△t TC 历元相对于滑动窗口起始的相对时间为t c,real ，则动基准在预报时间间隔△t p 对应的相对时间为t c,real +△t p ，此时动基准位置增量
Figure 92771DEST_PATH_IMAGE078
三维位置均通过下式计算得到：
Figure 1821DEST_PATH_IMAGE079
式中
Figure 353299DEST_PATH_IMAGE080
、
Figure 613379DEST_PATH_IMAGE081
...
Figure 957772DEST_PATH_IMAGE082
、
Figure 721329DEST_PATH_IMAGE083
、
Figure 758555DEST_PATH_IMAGE084
...
Figure 240352DEST_PATH_IMAGE085
和
Figure 404748DEST_PATH_IMAGE086
、
Figure 22812DEST_PATH_IMAGE087
...
Figure 230939DEST_PATH_IMAGE088
分别表示
Figure 200032DEST_PATH_IMAGE089
的x、 y和z三维方向的多项式系数，
Figure 886228DEST_PATH_IMAGE090
、
Figure 643619DEST_PATH_IMAGE091
和
Figure 22647DEST_PATH_IMAGE092
分别表示
Figure 213457DEST_PATH_IMAGE093
的x、 y和z三维方向在时间间隔△t p 内的动基准位置增量预报值。
10.根据权利要求1所述的卫星和惯性组合动对动实时精密相对定位方法，其特征在于：第五步中，合成相对定位结果dX RB,DGNSS/INS 为：
Figure 968924DEST_PATH_IMAGE094
。
CN202011327522.4A 2020-11-24 2020-11-24 一种卫星和惯性组合动对动实时精密相对定位方法 Active CN112147663B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011327522.4A CN112147663B (zh) 2020-11-24 2020-11-24 一种卫星和惯性组合动对动实时精密相对定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011327522.4A CN112147663B (zh) 2020-11-24 2020-11-24 一种卫星和惯性组合动对动实时精密相对定位方法

Publications (2)

Publication Number Publication Date
CN112147663A true CN112147663A (zh) 2020-12-29
CN112147663B CN112147663B (zh) 2021-02-09

Family

ID=73887426

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011327522.4A Active CN112147663B (zh) 2020-11-24 2020-11-24 一种卫星和惯性组合动对动实时精密相对定位方法

Country Status (1)

Country Link
CN (1) CN112147663B (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113009537A (zh) * 2021-02-18 2021-06-22 中国人民解放军国防科技大学 一种惯导辅助卫导相对定位单历元部分模糊度求解方法
CN113203414A (zh) * 2021-05-21 2021-08-03 北京交通大学 一种基于gps+bds ppp/imu紧组合的列车定位方法
CN113406670A (zh) * 2021-06-17 2021-09-17 哈尔滨工程大学 一种北斗非机动型广播星历故障实时监测方法
CN113703024A (zh) * 2021-07-16 2021-11-26 交通运输部水运科学研究所 一种基于北斗的小型船舶动态管理***及管理方法
CN115407371A (zh) * 2022-09-02 2022-11-29 中国人民解放军国防科技大学 基于PPP-B2b的实时高精度时间传递方法及装置
CN117168499A (zh) * 2023-09-04 2023-12-05 武汉大学 一种高频动态目标基准位置估测方法及计算机可读介质

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104536026A (zh) * 2015-01-08 2015-04-22 中国航空无线电电子研究所 一种动态对动态实时测量***
CN107894232A (zh) * 2017-09-29 2018-04-10 湖南航天机电设备与特种材料研究所 一种gnss/sins组合导航精确测速定位方法及***
CN108873034A (zh) * 2018-03-30 2018-11-23 广州海格通信集团股份有限公司 一种惯导辅助载波模糊度解算的实现方法
CN109917436A (zh) * 2019-04-28 2019-06-21 中国人民解放军国防科技大学 一种卫星/惯性组合实时精密相对动基准定位方法
US20190324154A1 (en) * 2018-04-24 2019-10-24 Novatel Inc. Global navigation satellite system (gnss) antenna data link
WO2019216261A1 (ja) * 2018-05-10 2019-11-14 ソニー株式会社 情報処理装置、情報処理方法およびプログラム
CN111578935A (zh) * 2020-05-08 2020-08-25 北京航空航天大学 一种利用惯导位置增量辅助gnss模糊度固定的方法
CN111965685A (zh) * 2020-07-07 2020-11-20 北京自动化控制设备研究所 一种基于多普勒信息的低轨卫星/惯性组合导航定位方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104536026A (zh) * 2015-01-08 2015-04-22 中国航空无线电电子研究所 一种动态对动态实时测量***
CN107894232A (zh) * 2017-09-29 2018-04-10 湖南航天机电设备与特种材料研究所 一种gnss/sins组合导航精确测速定位方法及***
CN108873034A (zh) * 2018-03-30 2018-11-23 广州海格通信集团股份有限公司 一种惯导辅助载波模糊度解算的实现方法
US20190324154A1 (en) * 2018-04-24 2019-10-24 Novatel Inc. Global navigation satellite system (gnss) antenna data link
WO2019216261A1 (ja) * 2018-05-10 2019-11-14 ソニー株式会社 情報処理装置、情報処理方法およびプログラム
CN109917436A (zh) * 2019-04-28 2019-06-21 中国人民解放军国防科技大学 一种卫星/惯性组合实时精密相对动基准定位方法
CN111578935A (zh) * 2020-05-08 2020-08-25 北京航空航天大学 一种利用惯导位置增量辅助gnss模糊度固定的方法
CN111965685A (zh) * 2020-07-07 2020-11-20 北京自动化控制设备研究所 一种基于多普勒信息的低轨卫星/惯性组合导航定位方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
YI DONG 等: "Tightly Coupled GNSS/INS Integration with Robust Sequential Kalman Filter for Accurate Vehicular Navigation", 《SENSORS》 *
王凌轩 等: "INS辅助的实时动态GNSS单频周跳探测", 《武汉大学学报》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113009537A (zh) * 2021-02-18 2021-06-22 中国人民解放军国防科技大学 一种惯导辅助卫导相对定位单历元部分模糊度求解方法
CN113009537B (zh) * 2021-02-18 2023-10-31 中国人民解放军国防科技大学 一种惯导辅助卫导相对定位单历元部分模糊度求解方法
CN113203414A (zh) * 2021-05-21 2021-08-03 北京交通大学 一种基于gps+bds ppp/imu紧组合的列车定位方法
CN113406670A (zh) * 2021-06-17 2021-09-17 哈尔滨工程大学 一种北斗非机动型广播星历故障实时监测方法
CN113703024A (zh) * 2021-07-16 2021-11-26 交通运输部水运科学研究所 一种基于北斗的小型船舶动态管理***及管理方法
CN115407371A (zh) * 2022-09-02 2022-11-29 中国人民解放军国防科技大学 基于PPP-B2b的实时高精度时间传递方法及装置
CN115407371B (zh) * 2022-09-02 2023-08-15 中国人民解放军国防科技大学 基于PPP-B2b的实时高精度时间传递方法及装置
CN117168499A (zh) * 2023-09-04 2023-12-05 武汉大学 一种高频动态目标基准位置估测方法及计算机可读介质
CN117168499B (zh) * 2023-09-04 2024-04-05 武汉大学 一种高频动态目标基准位置估测方法及计算机可读介质

Also Published As

Publication number Publication date
CN112147663B (zh) 2021-02-09

Similar Documents

Publication Publication Date Title
CN112147663B (zh) 一种卫星和惯性组合动对动实时精密相对定位方法
EP3805803A1 (en) Precise point position and real-time kinematic (ppp-rtk) positioning method and device
Xiaohong et al. Instantaneous re-initialization in real-time kinematic PPP with cycle slip fixing
AU2008260578B2 (en) Distance dependant error mitigation in real-time kinematic (RTK) positioning
CN111983654B (zh) 一种基于gnss的北极区域电离层相位闪烁因子构建方法
US9494693B2 (en) Method, apparatus, and system for determining a position of an object having a global navigation satellite system receiver by processing undifferenced data like carrier-phase measurements and external products like ionosphere data
CN107193028B (zh) 基于GNSS的Kalman相对定位方法
CN105758401A (zh) 一种基于多源信息融合的组合导航方法和设备
Li et al. Review of PPP–RTK: Achievements, challenges, and opportunities
CN111273687A (zh) 基于gnss观测量和机间测距的多无人机协同相对导航方法
AU2009330687A1 (en) Navigation receiver and method for combined use of a standard RTK system and a global carrier-phase differential positioning system
CN111965673A (zh) 基于多gnss的单频精密单点定位算法的时间频率传递方法
CN111427068B (zh) 一种动对动平台局域增强卫星a类星历故障完好性监测方法
CN115267855B (zh) 一种gnss-ins紧组合中异常值探测方法和平差定位方法
CN112526569B (zh) 一种惯导辅助卫导相对定位多历元逐级模糊度求解方法
Li et al. Predicting atmospheric delays for rapid ambiguity resolution in precise point positioning
EP2092362B1 (en) Phase based measurement corrections
Dong et al. Low-latency, high-rate, high-precision relative positioning with moving base in real time
Yoder et al. Low-Cost Inertial Aiding for Deep-Urban Tightly Coupled Multi-Antenna Precise GNSS
US6704650B1 (en) Technique for accurate distance and velocity calculations using the global positioning system (GPS)
CN115561796A (zh) 一种电网无人机巡检实时定位方法和***
Chen et al. A new cycle slip detection and repair method for single-frequency GNSS data
CN112630811A (zh) 一种实时ppp-rtk组合定位方法
Gunning Safety critical bounds for precise positioning for aviation and autonomy
CN115267858A (zh) 一种区域导航***辅助的精密单点定位方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant