RTKLIB 中的天线相位缠绕误差修正

之前一直没有细看精密定位（PPP）中的相位缠绕，不过我觉得还是应该至少从应用层面理解PPP中每项修正。因此抱着不求甚解的心情（没有仔细去研究右旋极化信号、偶极等概念），只从应用的角度来看看RTKLIB是如何进行天线相位缠绕修正的。

相位缠绕概念

根据各文献的引用来看，参考文献[1]是最经典的文献，但是没有找到下载，只能从其他博客和论文[2~6]得到对以下问题的一些理解。

• 相位缠绕是什么？
  发射天线和接收天线绕极化轴方向发生相对旋转时，测量的载波相位值会发生变化，这个现象称为相位缠绕。

• 相位缠绕对码的影响？
  相位缠绕只影响载波相位测量，不影响码相位测量，在非精密定位中可以忽略其影响。

• 相位缠绕对RTK的影响？
  相位缠绕的整周数部分可以被整周模糊度吸收，其不足整周的小数部分会造成误差。当天线只绕其极化轴旋转时，相位缠绕误差对各卫星的相位观测值影响均相同，在定位算法中和接收机相位钟差耦合而被吸收 。因此在静态短基线测量以及RTK测量中，通常忽略相位缠绕的改正（我理解应该是双差后就像钟差一样消掉了）。

• 相位缠绕对PPP的影响？
  对精密单点定位来说，影响不能忽略，因其影响可能会达到分米级（也有说厘米级别误差）。

相位缠绕误差计算方法

下面这张图中所有参数的解释可以参见论文[4]。我们所需要知道基础知识：GNSS卫星为了保证地面信号强度，其Z 轴指向地心，同时为了保证卫星能够得到更多的太阳能，X 轴处于地心、太阳所在平面且朝向太阳，Y 轴垂直与X、Z 轴构成右手坐标系。理解这一点，方便我们理解卫星本体坐标系的单位矢量是如何计算出来的。

根据参考文献[3]，相位缠绕的小数部分$\varphi$由下式计算：

式中：

• D’、D 分别为卫星和接收机天线的有效偶极矢量；
• K 为卫星至接收机的单位矢量；
• S^X、S^Y 分别为卫星本体坐标系中的 X 轴和Y 轴在地固系下的单位矢量；
• R^X 和 R^Y 为接收天线固定坐标系的北方向和西方向的单位矢量（参考[3]中说R^Y 是东向，但从RTKLIB源码以及参考[6]，应该是西向）；
• sign 为取符号运算符。

完整的相位缠绕改正Δϕ 是由其上个历元和本历元的小数部分共同组成。下式中，ROUND为四舍五入函数；Δϕ_prev 为前一个历元的相位缠绕改正值。

model_phw函数

• 所在文件：ppp.c
• 功能说明：根据模型计算天线相位缠绕误差
• 参数说明：

 args:  I   gtime_t 	time	:			GPS时
        I   int 		sat		:     		卫星
        I   const char *type	: 			antenna type（貌似没有使用）
        I   int 		opt		：          positioning options
        I   const double *rs	：         	卫星ECEF位置
        I   const double *rr	:          	接收机ECEF位置
        O   double 		*phw	:      		天线相位缠绕误差

• 函数源码：

/* phase windup model --------------------------------------------------------*/
static int model_phw(gtime_t time, int sat, const char *type, int opt,
                     const double *rs, const double *rr, double *phw)
{
    double exs[3],eys[3],ek[3],exr[3],eyr[3],eks[3],ekr[3],E[9];
    double dr[3],ds[3],drs[3],r[3],pos[3],cosp,ph;
    int i;
    
    if (opt<=0) return 1; /* no phase windup */
    
    /* satellite yaw attitude model */
    if (!sat_yaw(time,sat,type,opt,rs,exs,eys)) return 0;
    
    /* unit vector satellite to receiver */
    for (i=0;i<3;i++) r[i]=rr[i]-rs[i];
    if (!normv3(r,ek)) return 0;
    
    /* unit vectors of receiver antenna */
    ecef2pos(rr,pos);
    xyz2enu(pos,E);
    exr[0]= E[1]; exr[1]= E[4]; exr[2]= E[7]; /* x = north */
    eyr[0]=-E[0]; eyr[1]=-E[3]; eyr[2]=-E[6]; /* y = west  */
    
    /* phase windup effect */
    cross3(ek,eys,eks);
    cross3(ek,eyr,ekr);
    for (i=0;i<3;i++) {
        ds[i]=exs[i]-ek[i]*dot(ek,exs,3)-eks[i];
        dr[i]=exr[i]-ek[i]*dot(ek,exr,3)+ekr[i];
    }
    cosp=dot(ds,dr,3)/norm(ds,3)/norm(dr,3);
    if      (cosp<-1.0) cosp=-1.0;
    else if (cosp> 1.0) cosp= 1.0;
    ph=acos(cosp)/2.0/PI;
    cross3(ds,dr,drs);
    if (dot(ek,drs,3)<0.0) ph=-ph;
    
    *phw=ph+floor(*phw-ph+0.5); /* in cycle */
    return 1;
}

• 处理过程：

1. 首先调用sat-yaw函数，根据卫星的姿态模型计算出卫星本体坐标系X,Y方向的单位矢量exs、eys，即上面公式里的S^X、S^Y
2. 计算卫星至接收机的单位矢量；
3. 计算接收机天线在当地坐标系的北向、西向单位矢量；
4. 根据公式以及前一次的相位缠绕误差计算当前时刻相位缠绕误差

参考文献

[1] [Wu et al., 1993] Wu, J., Wu, S., Hajj, G., Bertiguer, W. and Lichten, S., 1993. Effects of Antenna Orientation on GPS Carrier Phase Measurements. Manuscripta Geodaetica. 18, pp. 91-98.
[2] GPS从入门到放弃（二十三） — 相位缠绕
[3] 王嘉琛, 刘根友, 郭爱智,等. 相对定位双差模型中的天线相位缠绕误差分析[J]. 武汉大学学报：信息科学版, 2020(1):9.
[4] 范曹明, 王胜利, 欧吉坤. GPS/BDS卫星姿态异常对PPP相位缠绕的影响及其改正模型[J]. 测绘学报, 2016, 45(10):7.
[5] https://gssc.esa.int/navipedia/index.php/Carrier_Phase_Wind-up_Effect
[6] Marreiros, João Paulo Ramalho. Kinematic GNSS Precise Point Positioning: Application to Marine Platforms. Diss. Universidade do Porto (Portugal), 2012.