GPS从入门到放弃（二十六） --- RTKLIB函数解析

GPS从入门到放弃（二十六） — RTKLIB函数解析

为了贴合这个系列的标题“从入门到放弃”，在入门之后现在就要放弃此方向了。虽然感觉遗憾，暂时也没有办法。在此附上此系列最后一篇，希望能给大家一些帮助。

此文中一些函数解析参考了 https://www.cnblogs.com/taqikema/p/8819798.html，在此表示感谢！

rtksvrthread

void *rtksvrthread(void *arg)

• 所在文件：rtksvr.c
• 功能说明：rtk服务线程。
• 参数说明：无
• 处理过程：

1. 检查svr->state，若为0，线程结束。
2. 从input stream 读取数据。
3. 写入到 log stream。
4. 若为精密星历，调用 decodefile 解码；否则调用 decoderaw 解码。
5. 对每一个观测数据，调用 rtkpos 进行定位计算。
6. 若定位成功，调整时间，写solution；若没成功，写solution。
7. 若有需要，发送 nmea request 给 base/nrtk input stream。
8. 休眠等下一个cycle。

rtkpos

int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)

• 所在文件：rtkpos.c
• 功能说明：根据观测数据和导航信息，计算接收机的位置、速度和钟差。
• 参数说明：

函数参数，4个：
rtk_t    *rtk      IO  rtk控制结构体
obsd_t   *obs      I   观测数据
int      n         I   观测数据的数量
nav_t    *nav      I   导航数据
返回类型:
int                O   (0:no solution,1:valid solution)

• 处理过程：

1. 设置基准站位置，记录观测值数量。
2. 调用 pntpos 进行接收机单点定位。若为单点定位模式，输出，返回。
3. 若为 PPP 模式，调用 pppos 进行精密单点定位，输出，返回。
4. 若无基准站观测数据，输出，返回。
5. 若为移动基站模式，调用 pntpos 进行基站单点定位，并加以时间同步；否则只计算一下差分时间。
6. 调用 relpos 进行相对基站的接收机定位，输出，返回。

pntpos

int pntpos (const obsd_t *obs, int n, const nav_t *nav, const prcopt_t *opt,
            sol_t *sol, double *azel, ssat_t *ssat, char *msg)

• 所在文件：pntpos.c
• 功能说明：依靠伪距和多普勒频移测量值来进行单点定位，给出接收机的位置、速度和钟差。
• 参数说明：

函数参数，8个：
obsd_t   *obs      I   观测数据
int      n         I   观测数据的数量
nav_t    *nav      I   导航数据
prcopt_t *opt      I   处理过程选项
sol_t    *sol      IO  solution
double   *azel     IO  方位角和俯仰角 (rad) (NULL: no output)
ssat_t   *ssat     IO  卫星状态            (NULL: no output)
char     *msg      O   错误消息
返回类型:
int                O   (1:ok,0:error)

• 处理过程：

1. 当处理选项 opt 中的模式不是单点模式时，电离层校正采用 broadcast 模型，即Klobuchar模型，对流层校正则采用 Saastamoinen 模型；相反，当其为单点模式时，对输入参数 opt 不做修改。
2. 调用 satposs 计算卫星们位置、速度、时钟
3. 调用 estpos 根据伪距估计接收机位置，其中会调用 valsol 进行 ${\chi }^2$ 检验和 GDOP 检验。
4. 若3中的检验没通过，调用 raim_fde 进行接收机自主完好性监测，判决定位结果的有效性，并进行错误排除。
5. 调用 estvel 根据多普勒频移测量值计算接收机的速度。
6. 赋值给卫星状态结构体ssat。

satposs

void satposs(gtime_t teph, const obsd_t *obs, int n, const nav_t *nav,
             int ephopt, double *rs, double *dts, double *var, int *svh)

• 所在文件：ephemeris.c
• 功能说明：按照所观测到的卫星顺序计算出每颗卫星的位置、速度、{钟差、频漂}
• 参数说明：

函数参数，9个：
gtime_t  teph      I   time to select ephemeris (gpst)
obsd_t   *obs      I   观测量数据
int      n         I   观测量数据的数量
nav_t    *nav      I   导航数据
int      ephopt    I   星历选项 (EPHOPT_???)
double   *rs       O   卫星位置和速度，长度为6*n，{
   x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      O   卫星钟差，长度为2*n， {
   bias,drift} (s|s/s)
double   *var      O   卫星位置和钟差的协方差 (m^2)
int      *svh      O   卫星健康标志 (-1:correction not available)
返回类型: 无

• 处理过程：

1. 大循环针对每一条观测数据，按顺序处理。
2. 首先初始化，将对当前观测数据的 rs、dts、var和svh数组的元素置 0。
3. 通过判断某一频率下信号的伪距是否为 0，来得到此时所用的频率个数。注意，频率个数不能大于 NFREQ（默认为 3）。
4. 用数据接收时刻减去伪距信号传播时间，得到卫星信号的发射时刻。
5. 调用 ephclk 函数，由广播星历计算出当前观测卫星与 GPS 时间的钟差 dt。注意，此时的钟差是没有考虑相对论效应和 TGD 的。
6. 用 4 中的信号发射时刻减去 5 中的钟差 dt，得到 GPS 时间下的卫星信号发射时刻。
7. 调用 satpos 函数，计算信号发射时刻卫星的位置(ecef,m)、速度(ecef,m/s)、钟差((s|s/s))。注意，这里计算出的钟差是考虑了相对论效应的了，只是还没有考虑 TGD。
8. 如果没有精密星历，则用广播星历的钟差替代。

• 注意：

1. 4中的公式原理：由 $\rho =c(t_u-t_s)$ 可得 $t_s=t_u-\rho /c$。

ephclk

int ephclk(gtime_t time, gtime_t teph, int sat, const nav_t *nav, double *dts)

• 所在文件：ephemeris.c
• 功能说明：通过广播星历来确定卫星钟差
• 参数说明：

函数参数，5个：
gtime_t  time      I   信号发射时刻
gtime_t  teph      I   用于选择星历的时刻 (gpst)
int      sat       I   卫星号 (1-MAXSAT)
nav_t    *nav      I   导航数据
double   *dts      O   卫星钟差，长度为2*n， {
   bias,drift} (s|s/s)
返回类型:
int                O     (1:ok,0:error)

• 处理过程：

1. 调用 satsys 函数，根据卫星编号确定该卫星所属的导航系统和该卫星在该系统中的 PRN编号。
2. 对于 GPS导航系统，调用 seleph 函数来选择最接近 teph 的那个星历。
3. 调用 eph2clk 函数，通过广播星历和信号发射时间计算出卫星钟差。

• 注意：

1. 此时计算出的卫星钟差是没有考虑相对论效应和 TGD的。

eph2clk

int eph2clk (gtime_t time, const eph_t *eph)

• 所在文件：ephemeris.c
• 功能说明：根据信号发射时间和广播星历，计算卫星钟差
• 参数说明：

函数参数，2个
gtime_t  time    I   time by satellite clock (gpst)
eph_t    *eph    I   broadcast ephemeris
返回类型：
double    satellite clock bias (s) without relativeity correction

• 处理过程：
  1. 计算与星历参考时间的偏差 dt = t-toc。
  2. 利用二项式校正计算出卫星钟差，从 dt中减去这部分，然后再进行一次上述操作，得到最终的 dt。(这一部分不知道是为什么？)
  3. 使用二项式校正得到最终的钟差。

satpos

int satpos(gtime_t time, gtime_t teph, int sat, int ephopt, const nav_t *nav,
           double *rs, double *dts, double *var, int *svh)

• 所在文件：ephemeris.c
• 功能说明：计算信号发射时刻卫星的 P(ecef,m)、V(ecef,m/s)、C((s|s/s))
• 参数说明：

函数参数，9个：
gtime_t time      I   time (gpst)
gtime_t teph      I   用于选择星历的时刻 (gpst)
int     sat       I   卫星号
nav_t   *nav      I   导航数据
int     ephopt    I   星历选项 (EPHOPT_???)
double  *rs       O   卫星位置和速度，长度为6*n，{
   x,y,z,vx,vy,vz}(ecef)(m,m/s)
double  *dts      O   卫星钟差，长度为2*n， {
   bias,drift} (s|s/s)
double  *var      O   卫星位置和钟差的协方差 (m^2)
int     *svh      O   卫星健康标志 (-1:correction not available)
返回类型:
int               O     (1:ok,0:error)

• 处理过程：

1. 根据不同的星历选项的值调用不同的处理函数，如果星历选项是 EPHOPT_BRDC，调用 ephpos 函数，根据广播星历计算出算信号发射时刻卫星的 P、V、C。如果星历选项是 EPHOPT_PREC，调用 pephpos 函数，根据精密星历和时钟计算信号发射时刻卫星的 P、V、C。

• 注意：

1. 此时计算出的卫星钟差考虑了相对论，还没有考虑 TGD。

ephpos

int ephpos(gtime_t time, gtime_t teph, int sat, const nav_t *nav,
           int iode, double *rs, double *dts, double *var, int *svh)

• 所在文件：ephemeris.c
• 功能说明：根据广播星历计算出算信号发射时刻卫星的 P、V、C
• 参数说明：

函数参数，9个：
gtime_t  time      I   transmission time by satellite clock
gtime_t  teph      I   time to select ephemeris (gpst)
int      sat       I   卫星号 (1-MAXSAT)
nav_t    *nav      I   导航数据
int      iode      I   星历数据期号
double   *rs       O   卫星位置和速度，长度为6*n，{
   x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      O   卫星钟差，长度为2*n， {
   bias,drift} (s|s/s)
double   *var      O   卫星位置和钟差的协方差 (m^2)
int      *svh      O   卫星健康标志 (-1:correction not available)
返回类型：
int                O    (1:ok,0:error)

• 处理过程:

1. 调用 satsys 函数，确定该卫星所属的导航系统。
2. 调用 seleph 函数来选择广播星历。
3. 根据选中的广播星历，调用 eph2pos 函数来计算信号发射时刻卫星的 位置、钟差和相应结果的误差。
4. 在信号发射时刻的基础上给定一个微小的时间间隔，再次计算新时刻的 P、V、C。与3结合，通过扰动法计算出卫星的速度和频漂。

• 注意：

1. 这里是使用扰动法计算卫星的速度和频漂，并没有使用那些位置和钟差公式对时间求导的结果。
2. 由于是调用的 eph2pos 函数，计算得到的钟差考虑了相对论效应，没有考虑 TGD。

eph2pos

void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts, double *var)

• 所在文件：ephemeris.c
• 功能说明：根据广播星历计算出算信号发射时刻卫星的位置和钟差
• 参数说明：

函数参数，5个
gtime_t  time      I   transmission time by satellite clock
eph_t    *eph      I   广播星历
double   *rs       O   卫星位置和速度，长度为6*n，{
   x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      O   卫星钟差，长度为2*n， {
   bias,drift} (s|s/s)
double   *var      O   卫星位置和钟差的协方差 (m^2)
返回类型：无

• 处理过程：

1. 通过卫星轨道半长轴 A 判断星历是否有效，无效则返回。
2. 计算规化时间 tk。
3. 根据不同卫星系统设置相应的地球引力常数 mu 和 地球自转角速度 omge。
4. 计算平近点角 M。
5. 用牛顿迭代法来计算偏近点角 E。参考 RTKLIB manual P145 (E.4.19)。
6. 计算升交点角距 u。
7. 计算摄动校正后的升交点角距 u、卫星矢径长度 r、轨道倾角 i。
8. 计算升交点赤经 O。
9. 计算卫星位置存入 rs 中。
10. 计算卫星钟差，此处考虑了相对论效应，没有考虑 TGD，也没有计算钟漂。
11. 用 URA 值来标定误差方差，具体对应关系可在 ICD-GPS-200H 20.3.3.3.1.3 SV Accuracy 中找到。

estpos

int estpos(const obsd_t *obs, int n, const double *rs, const double *dts,
           const double *vare, const int *svh, const nav_t *nav,
           const prcopt_t *opt, sol_t *sol, double *azel, int *vsat,
           double *resp, char *msg)

• 所在文件：pntpos.c
• 功能说明：通过伪距实现绝对定位，计算出接收机的位置和钟差，顺带返回实现定位后每颗卫星的{方位角、仰角}、定位时有效性、定位后伪距残差。
• 参数说明：

函数参数，13个：
obsd_t   *obs      I   观测量数据
int      n         I   观测量数据的数量
double   *rs       I   卫星位置和速度，长度为6*n，{
   x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      I   卫星钟差，长度为2*n， {
   bias,drift} (s|s/s)
double   *vare     I   卫星位置和钟差的协方差 (m^2)
int      *svh      I   卫星健康标志 (-1:correction not available)
nav_t    *nav      I   导航数据
prcopt_t *opt      I   处理过程选项
sol_t    *sol      IO  solution
double   *azel     IO  方位角和俯仰角 (rad)
int      *vsat     IO  卫星在定位时是否有效
double   *resp     IO  定位后伪距残差 (P-(r+c*dtr-c*dts+I+T))
char     *msg      O   错误消息
返回类型:
int                O   1表示成功，0表示出错

• 处理过程：

1. 初始化：将 sol->rr 的前 3 项位置信息（ECEF）赋值给 x 数组。
2. 开始迭代定位计算，首先调用 rescode 函数，计算当前迭代的伪距残差 $v$、几何矩阵 $H$、伪距残差的方差 var、所有观测卫星的方位角和仰角 azel、定位时有效性 vsat、定位后伪距残差 resp、参与定位的卫星个数 ns 和方程个数 nv。
3. 确定方程组中方程的个数要大于未知数的个数。
4. 以伪距残差的标准差的倒数作为权重，对 $H$ 和 $v$ 分别左乘权重对角阵，得到加权之后的 $H$ 和 $v$。
5. 调用 lsq 函数，根据 $dx=(H{H}^T{)}^{-1}$