@MinLily 2020-04-30T11:03:51.000000Z 字数 5452 阅读 3283

RTKlib相对定位源码解析：zdres函数

GNSS RTKlib 相对定位

最近阅读RTKlib开源代码，非常感谢“塔奇克敲代码”博主的博客（RTKLIB源码解析——单点定位），他将单点定位部分整理成函数小卡片，为我理解RTKlib提供了很大的帮助。我参照他的格式，记录整理对相对定位部分的个人理解。由于刚开始接触卫星定位，所以可能有理解不到位的地方，还请诸位指正。

zdres函数是rtklib相对定位relpos中的第二个函数。第一个函数satposs，“塔奇克敲代码”博主在单点定位中已有解释，因此不再赘述。而zdres函数中所调用的ecef2pos，satazel，satexclude函数，该博主同样也已进行阐述，因此也不再解释。下文暂未对潮汐修正函数tidedisp进行解释，但是RTKlib作者已经将参考论文列出，通过阅读论文应该也能自行理解。

我所基于的代码版本是RTKlib 2.4.3的一个拓展版本RTKexplore Demo5，这个版本主要针对低成本的GNSS定位。该版本整体算法并未做较大更改，只是针对低成本接收机进行了完善。

RTKlib相对定位源码解析：zdres函数
- - zdres
  - tropmodel
  - tropmapf
  - antmodel
  - zdres_sat

zdres

static int zdres(int base, const obsd_t *obs, int n, const double *rs,
                 const double *dts, const double *var, const int *svh,
                 const nav_t *nav, const double *rr, const prcopt_t *opt,
                 int index, double *y, double *e, double *azel)

所在文件：Rtkpos.c
功能说明：计算某历元所有卫星（未做差的）载波相位/伪距残差(残差 = 观测量 - 卫地距)
参数说明：

 args:  I   base:           1=base,0=rover 
        I   obs:             sat observations
        I   n:               # of sats
        I   rs [(0:2)+i*6]: sat position {x,y,z} (m)
        I   dts[(0:1)+i*2]: sat clock {bias,drift} (s|s/s)
        I   svh:            sat health flags
        I   nav:            sat nav data
        I   rr:             receiver pos (x,y,z)
        I   opt:            options
        I   index:          0=base,1=rover 
        O   y[(0:1)+i*2]:   zero diff residuals {phase,code} (m)
        O   e:              line of sight unit vectors to sats(卫星观测方向单位矢量)
        O   azel:           [az, el] to sats （方位角、仰角）                                    
return:
int                  O     (1:ok,0:error)

调用关系：

处理过程：
1. 将y数组中的所有元素初始化为0；
2. 如果没有接收机位置，返回0。如果是计算基站的残差，这里的位置输入参数是基站位置，如果计算移动站，位置输入参数是kalman滤波中的位置状态量；
3. 如果在配置中选择进行潮汐误差修正，调用tidedisp函数，对接收机位置进行地球潮汐误差项的修正；
4. 调用ecef2pos函数，将接收机位置从地球固定坐标系（WGS84）转化到大地测量坐标系（纬、经、高）；
5. 调用geodist函数，计算接收机到所有卫星的几何距离（卫地距）、单位观测矢量；然后调用satazel函数得到卫星方位角、仰角；
6. 通过satexclude函数判断该卫星是否需要被排除；
7. 利用输入参数dts，对3中获得的卫地距进行卫星钟差修正；
8. 调用tropmodel、tropmapf函数，对4中获得的卫地距进行流层延迟模型修正；
9. 调用antmodel函数，计算接收机天线偏移量；
10. 调用zdres_sat计算每颗卫星的载波相位/伪距残差，存放在y数组中。
注意事项：
1. zdres计算的相位/伪距残差，都是基于的原始观测量，并未进行做差；
2. 卫地距在整个过程中的修正包括潮汐修正（可选项）、卫星钟差修正、对流层修正，并计算了天线偏移误差。

tropmodel

extern double tropmodel(gtime_t time, const double *pos, const double *azel,
                        double humi)

所在文件：rtkcmn.c
功能说明：通过标准大气模型和saastamoinen模型计算对流层延迟（m）
参数说明：

* args   : gtime_t time     I   time
*          double *pos      I   receiver position {lat,lon,h} (rad,m)
*          double *azel     I   azimuth/elevation angle {az,el} (rad)
*          double humi      I   relative humidity
* return : tropospheric delay (m)

处理过程：
1. 先计算一些标准大气值，包括总压 $p$ 、大气温度 $T$ 、水汽压力 $e$ ，公式可见RTKlib manual E.5 (p149);
2. 对流层延迟包括了静力学延迟（即干延迟）和湿延迟两部分，该函数先计算了静力学延迟trph，然后计算了湿延迟trpw。
注意事项：
1. 该函数计算静力学延迟和湿延迟时所使用的公式和manual上不一致。但代码并没有错误，代码中使用的是正确的Saastamoinen模型公式：
  a) 可参考这篇博客。最原始的Saastamoinen的论文《Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites》由于找到的下载网站都需要收费，所以我参考了另一篇SCI论文《A Comprehensive Evaluation and Analysis of the Performance of Multiple Tropospheric Models in China Region》，也和代码中公式一致。
  
  $T_h = \frac{0.0022768p}{1.0 - 0.00266cos(2\phi)-2.8h\times10^{-7}}\times\frac{1}{cosz}$
  
  $T_w = 0.0022768(\frac{1255.0}T+0.05)e\times\frac{1}{cosz}$
  
  $T_r = T_h+T_w$
  b) 另一篇博客中提到《GPS测量与数据处理》第三版中是乘以(1.0-0.00266*cos(2.0*pos[0])-0.00028*hgt/1E3)，而不是除以。我感觉可能是书中的书写错误，大部分文献中还是用的除以。
  c) Saastamoinen模型公式貌似也有很多变形和近似，就没有再进行深入研究了。
2. 在zdres中使用该函数时，传入参数azel[]={0.0,90.0*D2R}，humi为0，因此实际上求取的是天顶方向的干延迟，接下来的tropmapf函数将计算天顶方向到观测方向的对流层延迟投影函数。

tropmapf

extern double tropmapf(gtime_t time, const double pos[], const double azel[],
                       double *mapfw)

所在文件：rtkcmn.c
功能说明：计算对流程延迟投影函数（即天顶方向到接收机相对卫星观测方向上的对流层延迟投影系数）
参数说明：

* args   : gtime_t t        I   time
*          double *pos      I   receiver position {lat,lon,h} (rad,m)
*          double *azel     I   azimuth/elevation angle {az,el} (rad)
*          double *mapfw    IO  wet mapping function (NULL: not output)
* return : dry mapping function

处理过程：
1. 这个函数中有两种投影函数的计算方法，分别是GMF和NMF，对应的两篇参考论文为“Global mapping functions for the atmosphere delay at radio wavelengths”和“Global Mapping Function(GMF): A new empirical mapping function base on numerical weather model data”。默认使用的是NMF方法，也可以通过定义IERS_MODEL宏来使用GMF方法。
2. 这个函数调用完成后，会将返回的干投影函数和tropmodel中计算得到的天顶方向对流层干延迟相乘，从而得到接收机相对卫星观测方向上的对流层延迟。
注意事项：
1. 干投影函数是通过返回值获得的，而湿投影是通过输入/输出参数mapfw获得的。在zdres函数中，本函数输入的mapfw=NULL，即不输出湿投影。
我的疑惑：
1. zdres函数调用这个函数时，为什么仅计算干延迟的投影函数？因为干延迟占延迟总量的90%，所以湿延迟忽略不计了吗？

antmodel

extern void antmodel(const pcv_t *pcv, const double *del, const double *azel,
                     int opt, double *dant)

所在文件：rtkcmn.c
功能说明：根据接收机天线的相位中心参数计算天线偏移量
参数说明：

 args   : pcv_t *pcv       I   antenna phase center parameters
          double *del      I   配置中所设置天线参考点（ARP）相对于地面标识的偏移
          double *azel     I   azimuth/elevation for receiver {az,el} (rad)
          int     opt      I   option (0:only offset,1:offset+pcv)
          double *dant     O   range offsets for each frequency (m)
 return : none
 notes  : current version does not support azimuth dependent terms

处理过程：
1. 生成一组接收机相对于卫星观测方向的单位矢量e，e[0]，e[1]，e[2]分别为东、北、天三个方向上的分量；
2. 频段不同，天线的相位中心偏移(PCO)不同。先计算出每个频段天线在东、北、天三个方向总的偏移，即相位中心偏移pcv->off[i][j]与del[j]之和；
3. 计算2中的相位中心偏移在观测单位矢量e上的投影dot(off,e,3)，由于e是单位矢量，所以和它求内积实际上就在该方向上的投影；
4. 计算天线相位中心变化量（PCV）：不同的高度角，相位中心变化不同，因此根据高度角对pcv->var[i]进行插值计算。
5. 3、4两部分即为总的天线偏移：dant[i]=-dot(off,e,3)+(opt?interpvar(90.0-azel[1]*R2D,pcv->var[i]):0.0);
注意事项：
1. 天线偏移量总共包括了三部分：
  a) 天线相位偏移（PCO） ARP –> APC（天线参考点到天线平均相位中心）
  b) 天线相位变化（PCV） APC –> IPC（平均相位中心到瞬时相位中心）
  c) 天线参考点（ARP）相对于地面标识的偏移
2. 相位中心偏移pcv->off[i][j]中的值来自于天线PCV文件。另外我在理解该函数的过程中，也看到另外两篇博客（博客1，博客2）对RTKlib中的天线偏移有所解释，如果不理解可以作为辅助阅读。

zdres_sat

static void zdres_sat(int base, double r, const obsd_t *obs, const nav_t *nav,
                      const double *azel, const double *dant,
                      const prcopt_t *opt, double *y)

所在文件：rtkpos.c
功能说明：计算单个卫星（未做差的）载波相位/伪距残差
参数说明：

 args:  I   base:           1=base,0=rover 
        I   r:              卫地距
        I   obs:            sat observations
        I   nav:            sat nav data
        I   azel:           [az, el] to sats （方位角、仰角）
        I   dant:           antenna offsets for each frequency (m)
        I   opt:            options
        O   y[(0:1)+i*2]:   zero diff residuals {phase,code} (m)

处理过程：
1. 如果在配置中选择了无电离层线性组合（IFLC：Ionospheric free linear combination）：
  a). 调用testsnr函数检查L1，L2观测量的载噪比是否大于配置中SNR Mask的最低要求（SNR Mask设置见RTKlib mannual 3.5章）；
  b). 计算无电离层线性组合系数C1、C2, 并采用该系数计算无电离层组合的天线偏移量；
  c). 计算无电离层线性组合载波相位和伪距残差：残差 = IFLC观测量 - 卫地距 - 天线偏移量。
2. 如果不是无电离层组合:
  a). 调用testsnr，检查每个频段的载噪比是否大于配置中SNR Mask的要求；
  b). 计算每个频段下的载波相位、伪距残差：残差 = 观测量 - 卫地距 - 天线偏移量。
注意事项：
1. 无电离层线性组合IFLC：
  
  $\rho_{1,2} = \frac{f_1^2 }{f_1^2-f_2^2}\rho_1 - \frac{f_2^2 }{f_1^2-f_2^2}\rho_2$