102-RTKLIB中的相位解缠

这是一篇未解决问题的博文

rtklib手册中的相位解缠公式：

首先$E_r$就很不明白，不知道为什么程序中要那样求解。
另外对于$E^s$，大小与前面卫星pco改正求得的各方向单位矢量大小相等，方向相反，猜测是pco中是坐标单位矢量，这里是转换矩阵，但是这里的$E^s$中xy的求解还是通过卫星的姿态进行求解的，具体没有细究，因为跟找到的资料对比还是没看明白。卫星的姿态相关计算可以参考下面几篇文章：
[1] 张宝成,欧吉坤,袁运斌,钟世明.GPS卫星姿态异常及其对精密单点定位估值的影响[J].科学通报,2010,55(Z2):2612-2618.
[2] 范曹明,王胜利,欧吉坤.GPS/BDS卫星姿态异常对PPP相位缠绕的影响及其改正模型[J].测绘学报,2016,45(10):1165-1170+1209.
[3] KOUBA J. A Simplified Yaw-attitude Model for Eclipsing GPS Satellites[J]. GPS Solutions, 2009, 13(1):1-12.

另外，$e_u^s$和$e_r^s$用的是同一个值，即ecef下接收机-卫星坐标的单位矢量，这个其实都一样，反正都是卫星到接收机天线方向的单位矢量。
大概就这些搞不透彻的。另外贴一下，

/* 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 */
	/* E.8.12星固系的参数 */
    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 */
	/* E.8.11中的参数 */
    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); /* eks=ek×eys 即E.8.13中 eus×eys */
    cross3(ek,eyr,ekr); /* ekr=ek×eyr 即E.8.14中 ers×eyy */
    for (i=0;i<3;i++) {
        ds[i]=exs[i]-ek[i]*dot(ek,exs,3)-eks[i]; /* E.8.13 */
        dr[i]=exr[i]-ek[i]*dot(ek,exr,3)+ekr[i]; /* E.8.14 */
    }
	/* E.8.15 没有加N的部分 */
    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;
}

/* satellite attitude model --------------------------------------------------*/
static int sat_yaw(gtime_t time, int sat, const char *type, int opt,
                   const double *rs, double *exs, double *eys)
{
    double rsun[3],ri[6],es[3],esun[3],n[3],p[3],en[3],ep[3],ex[3],E,beta,mu;
    double yaw,cosy,siny,erpv[5]={0};
    int i;
    
    sunmoonpos(gpst2utc(time),erpv,rsun,NULL,NULL);
    
    /* beta and orbit angle */
    matcpy(ri,rs,6,1);
    ri[3]-=OMGE*ri[1]; /* v1=vx-omega*ry */
    ri[4]+=OMGE*ri[0]; /* v2=vy+omega*rx , v3=vz */
    cross3(ri,ri+3,n); /* n=r×v */
    cross3(rsun,n,p); /* p=rsun×n */
    if (!normv3(rs,es)||!normv3(rsun,esun)||!normv3(n,en)||
        !normv3(p,ep)) return 0;
    beta=PI/2.0-acos(dot(esun,en,3));/* acos(-x)=pi-acos(x) */
    E=acos(dot(es,ep,3));
    mu=PI/2.0+(dot(es,esun,3)<=0?-E:E);
    if      (mu<-PI/2.0) mu+=2.0*PI;
    else if (mu>=PI/2.0) mu-=2.0*PI;
    
    /* yaw-angle of satellite */
    if (!yaw_angle(sat,type,opt,beta,mu,&yaw)) return 0;
    
    /* satellite fixed x,y-vector */
    cross3(en,es,ex);
    cosy=cos(yaw);
    siny=sin(yaw);
    for (i=0;i<3;i++) {
        exs[i]=-siny*en[i]+cosy*ex[i];
        eys[i]=-cosy*en[i]-siny*ex[i];
    }
    return 1;
}

对于程序中$\beta$的求解，程序中的计算方法为
$$\beta =\pi /2-arccos[(r_{sat}\times v)\cdot r_{sun}]$$

而相关文档中给出的公式是
$$\beta =arccos[(v\times r_{sat})\cdot r_{sun}]-\pi /2$$

可以证明其等价：
$$\begin{gathered} cos(A)=x \\ cos(\pi -A)=-x \end{gathered}$$

那么
$$\begin{gathered} arccos(x)=A \\ arccos(-x)=\pi -A \end{gathered}$$

即：
$$arccos(x)=\pi -arccos(-x)$$