1、 RTKlib 关于高精度 GPS 动态定位处理过程第一章引言 .31.1 调用主函数 main（ rnx2rtkp.c） 31.2 调用后处理函数 postpos（postpos.c） 31.3 处理基站信息 execses_b（postpos.c） .31.4 处理流动站信息 execses_r（postpos.c） .41.5 执行处理操作 execses（postpos.c） 41.6 函数调用流程图 .4第二章文件读取 .52.1 观测文件读取 readobsnav (postpos.c).52.1.1 文件头读取 redarnxh (rinex.c)52.1.2 文件的记录数据读取
2、 readrnxobs (rinex.c) 62.2 导航电文文件读取 .72.2.1 文件头读取 .72.2.2 文件的记录数据读取 .7第三章计算基准站位置和速度 .83.1 利用导航文件与基准站观测文件求卫星位置、速度和卫星钟钟差 satposs(ephemeris.c).83.1.1 卫星钟钟差计算 ephclk(ephemeris.c)83.1.2 卫星位置计算 satpos(ephemeris.c)93.2 码伪距单点定位 estpos(pntpos.c) 103.3 函数调用流程图 .11第四章动态相对定位求流动站位置 .124.1 码伪距单点定位求流动站的近似坐标 pntpos
3、 .124.2 载波相位动态相对定位 relpos(rtkpos.c) 124.2.1 利用导航文件和流动站观测文件求卫星位置和卫星钟钟差 satposs(ephemeris.c)124.2.2 求基准站对应的非差残差项 zdres(rtkpos.c).124.2.3 实时状态更新 udstate(rtkpos.c).134.2.4 求流动站对应的非差残差项 zdres(rtkpos.c).144.2.5 求双差残差项 ddres(rtkpos.c)144.2.6 卡尔曼滤波 filter(rtkcmn.c) 174.2.7 模糊度整数估计 resamb_LAMBDA().184.3 函数调用
4、流程图 20第五章总结 .215.1 结果输出 .215.2 不足之处 215.3 下一阶段计划与安排 21第一章引言精密 GPS 动态测量采用载波相位差分技术，其标准测量模式为，一台 GPS 接收机置于已知点，作为基准站来进行静态测量，另一台GPS 接收机置于载体上，作为流动站来进行动态测量。两台接收机同步观测相同的卫星，然后将两台 GPS 接收机的观测值进行组合处理，就可以获得流动站相对于基准站的坐标和速度。本文主要介绍用 RTKLIB 实现精密动态定位的过程，该过程包括观测文件和导航文件的读取、基准站位置的计算、流动站位置的求解、运行结果的输出，下面将分块讨论。1.1 调用主函数 mai
5、n（rnx2rtkp.c）读一个大型程序首先要找到主函数，rtklib 的主函数位于文件rnx2rtkp.c 中，主函数中主要定义了一些变量的缺省值。比如定位处理选项、解决方案输出选项等。1.2 调用后处理函数 postpos（postpos.c）输入文件包括观测文件、导航文件、精密星历文件等，postpos 在处理输入文件时有两种方法，一种是输入文件可以只包含关键词，然后通过函数处理，将关键词用时间、基准站编号、流动站编号等代替；另一种是直接调用输入文件的文件名，postpos 主要是来判断是哪一种输入方式，然后调用相应函数。1.3 处理基站信息 execses_b（postpos.c）精密
6、 GPS 动态测量，需 要利用基站信息进行差分定位，execses_b就是用来处理基站信息.1.4 处理流动站信息 execses_r（ postpos.c）execses_r 函数是用来判断调用输入文件时是调用的关键词还是直接使用了文件名，然后分别进行处理。1.5 执行处理操作 execses（postpos.c）从 execses 正式开始执行处理操作，这其中包括文件函数readobsnav，数据的处理函数 antpos、procpos 等。1.6 函数调用流程图注：方框里的字母表示调用的函数，方框里括号内字母表示调用函数所在文件。main（rnx2rtkp.c）readobsnav(po
7、stpos.c)execses（postpos.c）execses_r（postpos.c ）execses_b（ postpos.c）postpos（postpos.c）antpos(postpos.c)procpos(postpos.c)第二章文件读取2.1 观测文件读取 readobsnav (postpos.c)2.1.1 文件头读取 redarnxh (rinex.c)RINEX 格式的文件头用于存放与整个文件有关的全局性信息，位于每个文件的最前部，其最后一个记录为“END OF HEADER”. 文件头中，每一记录的第 6180 列为该行记录的标签，用于说明第 160列中所表示的内
8、容。文件头存放有文件的创建日期、单位名、测站名、天线信息、测站近似坐标、观测值数量及类型、观测历元间隔等信息。RTKLIB 中文件头读取的函数调用过程见图 2-1。图 2-1 文件头读取的函数调用过程readobsnav (postpos.c)redarnxh (rinex.c)redarnxfp (rinex.c)readrnxfile (rinex.c)readrnxt (rinex.c)2.1.2 文件的记录数据读取 readrnxobs (rinex.c)RINEX 格式的记录数据紧跟在文件头的后面，在观测文件中存放的是观测过程中每一观测历元所观测到的卫星及载波相位、伪距和多普勒等类型
9、的观测值数据等。在观测值文件中，所记录的载波相位数据的单位为周，伪距数据的单位为 m。观测值所对应的时标（即观测时刻）是依据接收机钟的读数所生成的，而不是标准的 GPS 时，因而在该时标中含有接收机的钟差。RTKLIB 中文件头读取的函数调用过程见图 2-2。图 2-2 文件记录数据的函数调用过程readobsnav(postpos.c)readrnxobs (rinex.c)readrnxfp (rinex.c)readrnxt (rinex.c)readrnxfile (rinex.c)2.2 导航电文文件读取2.2.1 文件头读取导航电文文件文件头的格式与观测文件头基本相同，但包含内容略
10、有不同。导航电文文件的文件头存放有文件创建日期、单位名及其他一些相关信息，另外，还有可能会包含电离层模型的参数以及说明 GPS 时与 UTC 间关系的参数和跳秒等。导航文件的文件头读取的函数调用过程与观测文件文件头读取流程相同，见图 2-1.2.2.2 文件的记录数据读取导航电文文件中存放的是所观测卫星的钟差改正模型参数及卫星轨道数据等。由于广播星历每 2h 更新一次，因此，在导航电文文件中可能会出现某颗卫星具有多个不同参考时刻钟差模型改正参数和轨道数据的情况。GPS 导航电文文件数据记录节中的内容为按卫星和参考时刻存放的各颗卫星的时钟和轨道数据，每颗卫星一个参考时刻的数据占8 行，第 1 行
11、为卫星的 PRN 号和该卫星钟的参考时刻及其改正模型参数，第 28 行为该卫星的广播轨道数据。导航文件的文件记录数据读取的函数调用过程与观测文件文件头读取流程基本相同，只是在调用最后一个函数时，导航电文文件调用的是函数 readrnxnav。第三章计算基准站位置和速度GPS 提供了三种观测量：伪距、载波相位和多普勒频移。本章主要讨论利用伪距单点定位的方法求取基准站接收机的位置。基准站位置计算调用的函数主要是 antpos 函数，然后再由 antpos 函数调用其他函数。用antpos函数求基站位置有三种方法，第一种是用伪距单点定位的方法求取，调用的是avepos函数；第二种是直接从位置文件中读
12、取，调用的是getstapos函数；第三种是直接从 rinex的文件头中读取。本文中采用的是第一种方法。由于观测文件是以历元的方式存取数据，所以在avepos 中要逐历元调用函数pntpos(pntpos.c)处理数据。Pntpos函数是处理单点定位的函数，该函数包含了卫星位置和钟差计算（satposs ） 、接收机位置和速度计算（estpos/estvel）等过程。3.1 利用导航文件与基准站观测文件求卫星位置、速度和卫星钟钟差 satposs(ephemeris.c)3.1.1 卫星钟钟差计算 ephclk(ephemeris.c)1、求取卫星信号发射时刻（3-1）()/j jitsatr
13、ecP式中， 表示第 颗星， 表示第 个测站，此处只有一个测站即基准jii站， 为 P1 码的伪距观测值， 为光速， 表示以接收机钟为时c()itrec间基准的卫星信号接收时刻， 表示以卫星钟为时间基准的卫星()jtsa信号发射时刻。2、用广播星历求取卫星钟钟差（3-2 ）()jtsatTOC（3-3）201)ftft（3-4 ）j式中， 为卫星钟的参考时刻， 分别为卫星钟偏差、卫TC012ff、 、星钟的漂移、卫星钟的漂移速度， 为卫星钟钟差，卫星钟钟差需s经多次迭代求得。3.1.2 卫星位置计算 satpos(ephemeris.c)1、计算观测瞬间卫星的平近点角 M（3-5）()jjjt
14、sat（3-6）jktoe（3-7 ）30()GMntkA式中， 为卫星钟的 GPS 时， 为星历的参考时刻， 为万有引st teGM力常数 与地球总质量 之积， 、 为广播星历给出的参数。Gn2、计算偏近点角 E（3-8 ）sin()1coeME解上述方程需用迭代法（也可用微分迭代法） 。3、计算升交距角 、卫星矢径 、卫星轨道倾角uri（3-9a）21sinarctoeEu（3-9b）()A（3-9c）0diitk（3-10a）sin(2)cos(2)uuC（3-10b）sin(2)cos(2)rrCuu（3-10c）in()cos(2)siiCu4、计算卫星在轨道面坐标系中的位置（3-1
15、1a）cosjxru（3-11b）injy5、计算卫星在 ECEF 坐标系中的位置 jjjXYZ、 、（3-12 ）0eLtktos（3-13a ）cscinjjjXxyL（3-13b）insjjjY（3-13c ）sjjZy注：卫星速度可由卫星位置差分求得。3.2 码伪距单点定位 estpos(pntpos.c)1、求伪距残差（即误差方程自由项 ）rescode(pntpos.c)l伪距原始观测方程为：（3-14 ）+c(-)djjjjjjitropinoPP式中， 为卫星到基准站的几何距222=jjjjiiiXYZ离， 为接收机钟差， 为卫星钟钟差， 为对流层延迟，ijjtrop为电离层延
16、迟， 为观测噪声。jiondjP将（3-14）线性化后可得：（3-15）0(k+lm)+c()djjjjj jjjjiiiitropinoPPXYZ:式中， ， 为坐标改正222000=jj jjii iiiXYZ:、 、项， ， ， 。0kjji0ljji0jji将含未知量的项移到等式右边，重组后有（3-0cd(k+lm)+cjjjjj jjjtropinoPiiiiPXYZ :16）式中含有四个未知量 。iiXYZ:、 、 、令伪距残差 ，可得误差方程式0cdjjjjjjjtropinoPlP（3-17）TlAV式中， 12nll12.01.0nkkllAm iiXYZc2、最小二乘法求方
17、程解 lsq(rtkcmn.c)（3-18 ）1()TXAPl式中， 为观测值权矩阵，可由对应的每颗星的观测误差求得。由（3-18 ）可求得基准站坐标 0ibiiiXXYYZZ3.3 函数调用流程图antpos(postpos.c)avepos(postpos.c)pntpos(pntpos.c)第四章动态相对定位求流动站位置在动态相对定位中，参考站 A 的接收机固定在基线向量已知的点上，流动站接收机是移动的，需要确定其在任意历元的位置。本章主要介绍相对定位的处理过程，调用的函数是 procpos(postpos.c.)。在函数 procpos 中，首先调用函数 inputobs（postpo
18、s.c ）函数来输入每一历元的观测数据和星历数据，然后再调用 rtkpos(rtkpos.c)函数来进行精密定位。下面主要讲述在函数 rtkpos 中处理数据的具体流程。4.1 码伪距单点定位求流动站的近似坐标 pntpos码伪距单点定位求流动站坐标的过程与第三章所述过程基本相同，只是观测文件改变，并且对于每一个历元需求一个流动站位置。4.2 载波相位动态相对定位 relpos(rtkpos.c)4.2.1 利用导航文件和流动站观测文件求卫星位置和卫星钟钟差 satposs(ephemeris.c)该过程与 3.1 中的卫星位置计算过程相同，只是改变了观测文件。4.2.2 求基准站对应的非差残
19、差项 zdres(rtkpos.c)（4-1）0,1(+d)1jjjj jbbtropbbycant（4-2）,2 2jjjj jtrsatpos estpos estvel（4-3）0,31(+d)1jjjj jbbtropbbyPcant（4-4）,422jjjj jtr式中， 、 分别表示载波 L1、L2 对应的观测量， 、 分别为1jbj 1P2码观测量， 、 对应 L1、L2 的天线相位改正。jdant2jt4.2.3 实时状态更新 udstate(rtkpos.c)1、实时位置更新 udpos(rtkpos.c)由 4.1 每个历元都会利用码伪距单点定位计算出一个流动站的近似位置，
20、再根据近似位置经过载波差分定位求出流动站的高精度坐标。2、周跳探测文中主要采用电离层残差法进行周跳探测，过程如下。组合观测量为：（4-5）21112jjjjj jIfffINcf上式消除了接收机至卫星的几何距离、轨道误差、接收机钟差、卫星钟差和对流层延迟误差，仅与电离层延迟、组合模糊的和观测值噪声有关。它与载体的运动状态无关，可用于静态或动态测量中零差、单差或双差载波相位观测值的周跳探测。在没有周跳时，它随时间变化缓慢。一旦有周跳产生，它就会产生比较显著的变化。因此其相邻历元的差值可以用来检测周跳，其差值为：（4-111()2()()2()jj j jjf fNtttt6）式中， 。12jjj
21、f3、码相位组合求模糊度的浮点解 udbias(rtkpos.c)GPS 伪距、载波相位和多普勒观测量表示如下：（4-7a ）1+c(-)djjjjjjitropinoPP（4-7b ）2jjjjjjitri（4-7c）11c(-)+jjjjjjji tropinoPNd（4-7d）22+jjjjjjji tri对（4-7a） 、 （4-7c ）分别进行站间差分可得：（4-8）1c(-)c(-)jjjjjjrbrbP+cjjrb11+NNjjjjjjjjjjrrbrrbbf f（4-9）11+jjjjjjbrr bff将（4-8） 、 （4-9 ）进行码相位组合可得：（4-10）N()()jj
22、jjjjrbrbrbP4.2.4 求流动站对应的非差残差项 zdres(rtkpos.c)（4-11a）0,1(+d)1jjjj jrrtroprycant（4-11b ）,2 2jjjj jrrrtrr（4-11c）0,3()jjjj jrrrtropryPc（4-11d）,4+djjjj jrrrtrrant式中， 、 分别表示载波 L1、L2 对应的观测量， 、 分别为1j2j 1P2码观测量， 、 对应 L1、L2 的天线相位改正。jdantjt4.2.5 求双差残差项 ddres(rtkpos.c)双差即不同测站，同步观测同一组卫星，所得的单差之间的差值。双差可以消除接收机钟差。另一
23、方面，双差放大了多路径效应和观测噪声的影响。1、双差方程（以 L1 为例）这里只考虑两颗卫星信号频率相同的情况，则双差方程的最终形式为：（4-12）,()()ijijijbrbrbrttNt式中， 。, ()ijjijirrrbt2、观测方程线性化（4-,0 ,()()()()Trijijijijijij rbrrbXYZjbirXYZtttattatN13）其中， ,0,0,0()()ijjirrrttt,ijjibb,iiirrN,jjjrbrb, ,0,0()()()+i jij rrXtXtttat, ,0,0()()()i jij rrYtYtt, ,0,0()()()+i jij
24、rrZtZtZat2223030k式中， 对应参考星， 对应参考站，参考站的坐标已知，ib为状态向量， 为状态向量的方差矩阵,jirrrrXYZNk（在程序中用 P 表示） 。因模糊度的浮点解已在 4.2.3 中求出，则可求得双差残差项（4-12a）1 11()()()()jiijjiijjrbrbrbrbvyyNN（4-12b ）2 222jiijjiijjrrrr（4-12c）3()(3)jiijjrbrbvyy（4-12d）4(4)jiijjrrby式中， 对应参考星， 对应参考站， 、 分别表示载波 L1、L2i 1jv2j对应的双差残差项， 、 分别表示码观测量 P1、P2 对应的双
25、差3jv4j残差项。实例中，设有 7 颗星，其中以 2 号星为参考星，以 L1 为例，可得出 6 个双差残差项，即 ，则 的方差1345611vvv1协方差矩阵为（4-13）6000111555RijiRijiiij 式中， ，22()sinjbRjkacdel22()sinibRikacdel对应参考星， 对应卫星的高度角， 为载波相位观测误差，iel ba、为基线向量观测误差， 为卫星钟稳定误差。cd系数矩阵 和 状态向量 分别为：HX,1,2,6, ,106()()iiiXYYiiiZZattatttt1,2,310,4,5,6,7,rrbrbrbrbXYZNXN2210 23030k式
26、中， =2 表示 2 号星为参考星， 为状态向量 的方差矩阵（在i kX程序中用 表示） 。P4.2.6 卡尔曼滤波 filter(rtkcmn.c)1、模型方程为：（4-14）,1kkkXwlAe式中， 为系统误差向量， 为观测误差向量， ，k ke(0,)kkwN:。(0,)kkeNR:实例中，令 ， ， ， ， 为单位阵。,1kXTkHA1kvekR,1k则有：（4-15）TkQHR（4-16 ）1K（4-17）1PXKv（4-18）()TkIH式中， 对应浮点解。PX4.2.7 模糊度整数估计 resamb_LAMBDA()整周模糊度的固定可以提高基线向量解的精度。在载波相位数据处理中
27、，使用双差而非单差非常重要，因为在双差观测方程中，钟差参数已被消除，因此模糊度分离成为可能。整周模糊度成功解算的另一个重要因素是卫星的几何构型，在这里不做讨论。整周模糊度的解算主要有三个步骤。第一步是从算法角度考虑生成所有可能的整周模糊度组合；第二步是搜索出正确的模糊度组合；第三步是模糊度整数解的检验。下面主要讲述最小二乘模糊度降相关平差法即 lambda 算法求解整周模糊度。1、对 4.2.6 中所得结果进行变换得（4-19）PyDX（4-20）yQ式中，10911.0.1.0.1.0D 12,3,914,5,16,7,rrbrrbrrbrrXYZNyN（4-21）12,34,1566,7r
28、brrrbrrNN31rrXYZ为状态向量 的协方差矩阵，再对 做变换得 、 ， 对应yQyyQbaQ状态向量中模糊度浮点解， 为状态向量中模糊度的协方差矩阵，b为状态向量中非模糊度参数与模糊度的协方差矩阵。ab2、求模糊度的整数解对 进行 Cholosky 分解，即bQ（4-20 ）TNLD式中， ， ， 是单位下三角矩阵, 但非对角元不是整数，bNL对 取整为 ，做变换1（4-20）(1)1TNQ再对 做 Cholosky 分解，将分解的单位下三角矩阵取整，用 表(1)NQ 2L示，再做变换 （4-21）(2)(1)2TNL如此循环，直至分解的单位下三角矩阵取整后为单位阵，得到（4-22）
29、11TkkZNQLL 取 （4-23）k令 ，则 ，从而 ，得出模糊度的固定解 。1T Z1NTZ bN3、计算整周模糊度条件下的非模糊度参数解将整周模糊度 回代，计算模糊度条件下的非模糊度参数解 ：b X（4-24）131()aXQ（4-25）Tba式中， 对应状态向量中模糊度浮点解， 为状态向量中模糊度间bQ的协方差矩阵， 为状态向量中非模糊度参数与模糊度的协方差矩abQ阵，为非模糊度参数的浮点解。TrrXYZ4、修复单差模糊度由（4-21 ）可知，向量 表示的是双差模糊度解，将已解算出的整周b模糊度 回代b12,34,156,7rbrrrbrrNN若已知参考星所对应的单差模糊度即 ，便可
30、计算出每颗星对应的1,rb单差模糊度解。4.3 函数调用流程图第五章总结5.1 结果输出结果输出分为两部分，第一部分是输出文件头，在函数 execses 中调用 outhead(postpos.c)函数；第二部分是输出数据的处理结果，在函数 procpos(postpos.c)中调用函数 outsol(solution.c)。输出结果依次为历元所对应的时间、坐标、坐标方差与协方差、差分龄期、模糊度测试比率。procpos(postpos.c)rtkpos(rtkpos.c)relpos(rtkpos.c)satposs(ephris.c)zdres(rtkpos.c)udstate(rtkpos.c)zdres(rtkpos.c)ddres(rtkpos.c)filter(rtkpos.c)resamb_LAMBDA5.2 不足之处文中只是讲述了 RTKlib 用于动态定位的整个流程，许多细节问题未予深入探讨，比如周跳探测的方法、电离层延迟改正、对流层延迟改正均未详细讨论，双频载波观测值的线性组合也未涉及，整周模糊度的解算也只涉及到 LAMBDA 一种算法，要做的工作还很多，要学的东西还很多。5.3 下一阶段计划与安排下一阶段，会在前面工作的基础上，改进不足，对于周跳探测、模糊度解算的方法研究上多下功夫，同时在 RTKlib 与北斗的结合应用上多下功夫。
