帮“不想偏航”转发的专栏《RTKLIB代码学习》

记录一下这俩礼拜从零开始学习RTKLIB的一些经验，据说RTKLIB是功能好用，代码难读。如能帮助新入门的同道最好不过。如果哪里不对还请大神不吝在评论区指正。

调试代码使用工具为VS2019，RTKLIB版本为2.4.3 b33，是2018.1.30最后更新的源代码。

下载方面，国内从GitHub上下载比较慢，码云(Gitee)上有源码和软件包，下载比较快。

一、编译运行

RTKLIB的编译基本参照这个CSDN上的分享完成的https://blog.csdn.net/sd28you28/article/details/82911273

通常来说，大多数国内学习者都是将rnx2rtkp.c作为主函数，在此基础上进行修改和运行的。路径为..\RTKLIB\app\rnx2rtkp\rnx2rtkp.c        在这个路径里，msc文件夹里有一个msc.sln，这就是VS的解决方案，我们几乎可以直接使用了。当然，稳妥起见，最好把整个rtklib都备份一下，以免改着改着就崩了还忘了怎么改回去。打开后，如果你电脑的VS版本更高，通常会弹出升级工具包之类的窗口，按指示升级即可。以下走的是另外一条路，可供参考。

新建一个项目去做也可，我个人是新建了项目命名为PPPtest，然后拉进来了如下图所示的文件：

新建项目PPPtest引入的文件

然后选择Debug模式，在右边的解决方案资源管理器中，源文件下找到rnx2rtkp.c，双击点开，我是按照https://blog.csdn.net/sd28you28/article/details/82911273的指示把程序调到可以运行的。感谢原作者，在此不再赘述。

rnx2rtkp.c可作修改，下图是我试图作精密单点定位（precise point position）时做的一些修改：新增了字符串指针optfile指向由rtkpost.exe生成的conf文件，argc设置为1跳过了下面很多的循环，然后选择了定位模式(prcopt.mode)、卫星系统(prcopt.navsys)，频点选择(prcopt.nf ,3为L1+L2+L5)，导入了观测数据(.obs)、广播星历(.nav)、精密星历(.sp3)，精密钟差(.clk)，注意要ret=postpos...前面的4行if判定，因为之前的一大串把n赋为0了，留着最后的if会直接return -2，后面就没法运行了。

rnx2rtkp.c

一些关于PPP的修改

在运行中，由于我的星历文件里只有北斗的数据，运行时总是说no nav data就退出了，通过逐步调试发现，问题在MINPRNCMP和MINPRNMAX都赋了0（rtklib.h），在【项目->属性->C/C++->预处理器->预处理器定义->编辑】里，加上【ENACMP】（跟之前链接里处理格洛纳斯的报错的方式一样），就启用了北斗的相关数据。然后就可以运行了。

但是处理数据时发现，同样的观测数据、星历、设置，使用源代码运行的结果和使用rtkpost.exe运行的结果有差异，可见下图：

源代码运行结果和rtkpost.exe结果有差异

查看修改时间可以知道，rtkpost.exe更新时间是2019.8.19，而源码在2018年初就不再更新了，中间做了什么修改不得而知，亦或是更新的代码放在别处了？

二、精密星历读取

精密星历是区别于广播星历的星历文件，许多研究中心都通过获取世界各地的测站的数据，分析计算推出了各自的精密星历/钟差产品。就我所下载少量个精密星历而言，通常对于GPS的精密星历，各类产品的差异都比较小，定位误差大约在厘米级到分米级，钟差差异在几个纳秒(ns)左右。而对于北斗，差异是比较大的，对于北斗的MEO卫星，精密星历的定位误差通常都比较小，与GPS类似，但是IGSO就要差一些，GEO差异最大，在米级左右，相隔36000公里，一动不动的定位着实不易。钟差的差异通常存在系统性误差，来源大概是各机构的时间基准等不一致。这个系统误差精密单点定位中是会被接收机钟差吸收掉的，因此影响有限。

精密星历读取的调用流程为：main里的postpos->(openses中读取一些修正参数，如天线相位中心修正等)->execses_b->readpreceph->readsp3，在readsp3中读取精密星历和钟差，赋值到一个很大的结构体nav中，这个nav指向的是一个静态全局结构体navs，阅读它的定义会发现。。简直包罗万象，让我感觉之前自己写的程序都是幼儿园水平。

readsp3中，readsp3h函数读取文件头，第一行str2time把日期转换为它的时间计数方式，rtklib的时间计数既不是GPST也不是utc，而是把1970年1月1日0时0分0秒作为起点，不闰秒不跳秒。

文件preceph.c中的readsp3h函数

82行处有所修改，原来是ns=(int)str2num(buff,4,2)，这是因为近年导航卫星发射频繁，原来只要两位数就够用，现在已经达到一百多颗卫星，需要三位数存储了。事实上，访问北斗的一些官方网站，比如这个【http://www.csno-tarc.cn/system/constellation】。

84-88行是读取精密星历播发的卫星编号，每3个字符一颗卫星，94-95是读取星历中位置和钟差的标准差参数，进行误差估计与预报。但是这是为只有GPS的sp3文件设计的，不同的sp3文件里，这个参数不一定在第15行，这时就读不了了。

（这里本该有readsp3b函数截图，但是太长我就不放了）

readsp3b函数读取文件体，EOF是文件结尾字符，读取即跳出。*开头的一行都是时间，因此调用str2time转换为软件内通用的时间制式。141-152行为每个时间的星历生成存储空间并初始化为0。每个时间的星历前4为表示编号，第一位为P表示后面的是位置和钟差，据说，若第一位为V则表示后面的是速度和钟差变率，不过V开头的精密星历我还没下载到过。157-162行是把带字母的编号转换成数字顺序，毕竟数组调用时不认字母。接下来便是对位置和钟差信息验证和赋值，在199行调用addpeph函数把读取完的一个时刻的精密星历结构体peph赋给nav结构体。

在execses_b中，readpreceph接下来是读取基站名称之类的操作，个人感觉不是很重要，然后是读取扩展路径，进入execses_r函数，r表示rover，移动站，与之前的b（base，基站）区分，然后也是读取流动站名称之类的操作，然后进入execses函数，这个函数囊括了读写和运算的顺序，本文在第四节将解释其中的setpcv部分。在观测数据、广播星历、天线相位修正、海洋潮汐负荷读取结束，调试日志（.trace）和定位结果文件（.pos）文件头写入之后，进入处理模式选择，包括前向、后向、联合三种，和rtkpost.exe在options中第一页的选择一致。最终都会进入procpos函数进行处理。在这个函数中完成模式和卫星系统的一些选择，又指向rtkpos函数(文件rtkpos.c中)进行处理。在这里，根据之前设置的选择配置，开始进入不同的定位函数，相对定位的、单点定位的等等，这里假设进入pppos函数（ppp.c文件）。卫星位置和钟差计算在satposs函数中(ephemeris.c)。

satposs函数里，先读取了当前接收机观测时刻里所有的伪距信息，并由此除以光速，加上广播星历中的钟差改正，计算出了卫星信号播发的时刻。由于有精密星历，进入satpos函数后，在swtich选择中会进入peph2pos函数，这就是计算给定时刻（即卫星信号播发时刻）的由精密星历计算的卫星位置和钟差的函数。

三、精密星历下卫星位置和钟差的计算

在peph2pos函数（preceph.c文件）中，相关参数的计算流程如下：先计算卫星信号播发时刻的位置和钟差（调用了pephpos，如有clk文件还会调用pephclk文件），然后计算0.001秒后的位置和钟差，两个位置作差除以时间间隔即是速度（求导麻烦，编程中通常都靠自变量增加一个小量然后作差除以小量来近似计算导数），然后进行卫星天线相位中心的计算（精密星历确定的是卫星几何中心的位置，但是信号播发是天线完成的，这中间有一个杆臂矢量需要修正），如有精密钟差产品还会进行相对论效应的修正。

接下来学习pephpos函数，satantoff将在第四节叙述。

pephpos函数的一部分

preceph.c文件405-410行，是对比了卫星信号播发时刻是否在读取的精密星历起止时间范围内。因为rtklib里对星历的插值是用的多项式插值法，数值分析中我们知道，高次插值会发生龙格现象，延伸段的结果可信度较低。因此如果给定的信号播发时刻距离精密星历起止太远，会显示没有精密星历可用。

412-415行用的是二分法搜索，搜索距离卫星信号播发时刻最近的星历数组指针。三行代码就凝练地完成了二分法搜索功能，再次告诉我我就是个辣鸡，还是不可回收的那种。

419-428行算是在校验参数，因为读取的是结构体中的数据，没有读取到的或者无效的部分都赋值为0了，，因此这一部分的观测值时没有精密星历数据的，需要跳过。

pephpos函数的一部分

429-442行进行了地球自转改正，卫星位置是使用的ECEF系进行表示，因此z轴无需改正。

443-445行是调用了interppol函数进行多项式插值计算，使用了Neville算法，默认是11个节点进行10阶多项式插值。其实我也不是很懂为什么要用这么高的次数，但是插值精度确实很高，尤其是满足前后超过5个节点时。

446-454行是对位置的计算结果进行误差估计。

455-459行是对钟差线性插值参数的计算，在472行对信号播发时刻的钟差进行了线性插值。

pephclk函数对于钟差的计算也是如此，但是clk文件钟差给的比较频繁，比如30s左右就给一个，比起sp3文件5分钟、15分钟给一个要频繁得多，而钟差通常变化比较慢，因此另外的clk文件确实能够提供更高的钟差插值精度。

四、卫星天线相位中心修正

rtklib手册对于卫星天线相位修正和体坐标系的示意图

卫星天线相位修正通常是三步：一是读取atx文件，二是计算天线相位中心偏移（pco）的修正量，三是根据地面接收机的位置计算相对于卫星天线相位中心的方位，进行天线相位变化（pcv）的修正。由于atx文件中没有北斗卫星的pcv信息（都是0），因此本节只作pco修正的叙述，敬请谅解。

卫星天线相位中心的读取和设置的调用路径为：main里的postpos->openses->readpcv，atx文件到现在十几兆，很多并不是卫星的相位中心修正，而是接收机的修正。接收机的修正也是调用了这个readpcv函数。对于atx文件，readpcv里通过readantex函数（rtkcmn.c文件中）进行读取，存储在pcvs结构体中，这个结构体最终指向的也是一个静态全局结构体pcvss，相应的，接收机也有个pcvsr结构体存储信息。

在atx文件中，第60个字符后有该行的说明（如果该行是一大串数据就没有）。START OF ANTENNA表示开始显示一次新的天线相位中心修正信息，到END OF  ANTENNA。由于卫星存在更替、相位测量结果存在变化，因此同一颗卫星可能有好几个相位中心修正信息，这个由有效期来控制，VALID FROM和VALID UNTIL分别给出了相位修正的有效起止时间。各个频点的相位中心可能不一致，由START OF FREQUENCY和END OF FREQUENCY控制。这里也有个我比较奇怪的地方，就是下面截图的C01星，是北斗二代的地球静止卫星，播发三个频点的信号（B1I/B2I/B3I），怎么会有4个频点信息呢？如有大神路过还请不吝指教。NORTH/EAST/UP是频点的天线相位中心偏移，NOAZI是天线相位中心变化，隔一定角度给出一个值，这个角度由之前的ZEN1/ZEN2/DZEN计算给定。

igs14_2118.atx文件部分内容截图

readantex函数的源代码就不全放了，主要是根据atx文件60位之后的说明信息不同，将各种信息存储至pcv结构体中，还是比较好懂的。比较搞不懂的一点在2089-2092行，见下图。一颗卫星的载波信号总数NFREQ为3，但是通过sscanf函数扫描近来的f只比较了前三个freqs（也即1,2,5号频点），那么北斗的频点在2,7,6号，就只能读取一个2号频点的（还有个不知怎么就有的1号频点的），那么这样的修正值显然是不对的，这让我百思不得其解，还请指教。

readantex函数部分截图

另外，访问【http://www.csno-tarc.cn/support/satelliteparameters】可以获得更加准确的卫星天线相位中心的偏移信息，不知道为什么atx文件迄今仍不采纳。插句题外话，据说国外的北斗使用效果不如GPS的原因，一方面是卫星更少，另一方面就国外设备不上心，从硬件制造时信号捕获的能力不足，到软件编写的方面。

卫星pcv信息在以年为单位时时频繁变化的，但在实验时就那么几天甚至几分钟，在这段时间通常是没有变化的。因此给定一个时间后就可以设置当前需要的pcv信息，并读取到navs结构体中。这里用到的函数就是setpcv，调用路径为main里的postpos->execses_b->execses_r->execses->setpcv。这个函数在postpos.c文件中，其中第826-834行是卫星天线pcv信息的设置。调用searchpcv函数读取pcv，然后放置进nav结构体中。

对于searchpcv函数，也是比较容易读懂的，不放源代码截图了。主要就是通过卫星prn号和修正信息的起止时间，逐个比较，确定需要pcv的信息并返回。

最后是对于读取的信息的使用，也即天线相位中心偏移的修正。读取atx文件所得的数据是卫星体坐标系下的表示，为了计算需要必须投影至ECEF系下表示，这就用到了satantoff函数（preceph.c文件中），它在获得了精密星历的位置和钟差结果之后进行修正，调用顺序在第三节中有。

satantoff函数截图

初看这个函数我是懵逼的，满脸这啥玩意儿啊咋回事啊，咋还和太阳位置扯上关系了呢。回到上面rtklib手册对pco信息的示意图，下面就是坐标系的解释，卫星发送信号的天线始终对地定向，据说通过调姿始终保持对地心±1°的偏差内，而同时，大功率信号发送需要足够的电能，因此绕卫星质心到地心的连线旋转，需使太阳能帆板尽量正对太阳。卫星对于活动部件其实是极为抗拒的，尤其是这种要服役十几年的卫星，太阳能帆板是不会活动的，只能依靠调姿进行，因此我们就可以假设，卫星体坐标系一轴指向地心，一轴指向太阳。这就是sunmoonpos函数的作用了，给出了ECEF系下卫星到太阳的方向矢量，和卫星到地心的矢量叉乘，即得到了体坐标系在ECEF下的表示，最后加上电离层修正（这一块原理我不大懂），就得到了卫星天线相位中心偏移量在ECEF下的表示，用dant表示。

接下来康康sunmoonpos函数，时间要转换为utc（感觉过不久就应该调整差值为19秒了呢），再转为ut1，然后计算太阳和月亮方向在地心惯性系下的表示，调用了sunmoonpos_eci函数，这一步触及了知识盲区。。各种计算我一个都看不懂。跳过这一步往下看看，是计算地心惯性系（ECI）到地心地固系（ECEF）的转换函数，然后坐标系转换矩阵乘矢量，就得到了新坐标系下的矢量表示。

从头懵到尾的eci2ecef函数

那么是如何转换的呢？打开eci2ecef函数看一下。。先是世界时、GPS时间系统和格林尼治平恒星时的混乱关系。。然后计算了个不知道干啥用的天文参数（不负责任的大胆猜测可能是极移的参数），然后是使用了iau1976岁差模型，又使用了iau1980章动模型，计算了格林尼治真恒星时，最后就一通操作得到了该时刻ECI到ECEF系的转换矩阵和格林尼治平恒星时（假装看懂了）。至于中间的岁差、章动、极移（猜测）的函数，有兴趣可以继续看看，我是已经看不懂了。

好了，以上就是这篇学习经历的全部内容了。不得不承认自己的水平实在有限，文中谬误之处还请各位看官指正。rtklib作为开源软件到现在已经完善了十几年了，感谢全世界参与其中工作的学者们，身为后辈希望能够尽早追上前人的脚步，在学术上做一点微不足道的贡献。