rtklib-eph2pos-利用广播星历计算卫星的PVC-详细解说

---

前言

rtklib–伪距单点定位（single-point positioning）学习（1）Satposs这篇文章中介绍了satposs整体的流程，但是当时因为篇幅限制原因，没有对如何利用广播星历中的轨道参数对卫星进行定轨进行介绍，只说了大致的流程，也就是eph2pos这个函数没有展开介绍。很多同学应该会遇到给卫星定轨的编程问题，所以这里详细介绍一下，之后会介绍利用精密星历进行卫星定轨的编程实现。

/* broadcast ephemeris to satellite position and clock bias --------------------
* compute satellite position and clock bias with broadcast ephemeris (gps,
* galileo, qzss)
* args   : gtime_t time     I   time (gpst)
*          eph_t *eph       I   broadcast ephemeris
*          double *rs       O   satellite position (ecef) {x,y,z} (m)
*          double *dts      O   satellite clock bias (s)
*          double *var      O   satellite position and clock variance (m^2)
* return : none
* notes  : see ref [1],[7],[8]
*          satellite clock includes relativity correction without code bias
*          (tgd or bgd)
*-----------------------------------------------------------------------------*/
extern void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts,
                    double *var)
{
    double tk,M,E,Ek,sinE,cosE,u,r,i,O,sin2u,cos2u,x,y,sinO,cosO,cosi,mu,omge;
    double xg,yg,zg,sino,coso;
    int n,sys,prn;
    
    trace(4,"eph2pos : time=%s sat=%2d\n",time_str(time,3),eph->sat);
    
    if (eph->A<=0.0) {
        rs[0]=rs[1]=rs[2]=*dts=*var=0.0;
        return;
    }
    tk=timediff(time,eph->toe);

       首先判断半长轴是否小于0，如果是，那么卫星位置、钟差、钟漂、方差均设置为0，，然后tk算出信号发射时刻与星历参考时间的时间差。

计算卫星运动的平均角速度n与平近点角M

switch ((sys=satsys(eph->sat,&prn))) {
        case SYS_GAL: mu=MU_GAL; omge=OMGE_GAL; break;
        case SYS_CMP: mu=MU_CMP; omge=OMGE_CMP; break;
        default:      mu=MU_GPS; omge=OMGE;     break;
    }
    M=eph->M0+(sqrt(mu/(eph->A*eph->A*eph->A))+eph->deln)*tk;

调用satsys函数返回卫星的种类与prn号，如果是GPS卫星那么引力常数mu=GM=3.9860050E14，地球自转角速度omge=7.2921151467E-5 rad/s,sqrt(mu/(eph->Aeph->Aeph->A))计算出卫星的平均角速度，加上eph->deln摄动参数得出计算时刻的平均角速度。然后乘以tk得到观测时刻与参考时刻的平近点角的变化值，加上参考时刻toe的平近点角得到观测时刻的平近点角。

计算卫星运动的偏近点角E

for (n=0,E=M,Ek=0.0;fabs(E-Ek)>RTOL_KEPLER&&n<MAX_ITER_KEPLER;n++) {
        Ek=E; E-=(E-eph->e*sin(E)-M)/(1.0-eph->e*cos(E));
    }

这里使用了牛顿迭代法，求解 E=M+esinE 这个方程。

计算卫星运动的真近点角f与升角角距

 sinE=sin(E); cosE=cos(E);
 u=atan2(sqrt(1.0-eph->e*eph->e)*sinE,cosE-eph->e)+eph->omg;

u=atan2(sqrt(1.0-eph->e*eph->e)*sinE,cosE-eph->e)这是计算出卫星运动的真近点角，然后加上omg

计算卫星初始升交角距、矢径、轨道倾角

u=atan2(sqrt(1.0-eph->e*eph->e)*sinE,cosE-eph->e)+eph->omg;
r=eph->A*(1.0-eph->e*cosE);
i=eph->i0+eph->idot*tk;

这一步相对比较简单就是利用基本定义计算出一些基本参数。

对卫星初始升交角距、矢径、轨道倾角进行摄动改正

sin2u=sin(2.0*u); cos2u=cos(2.0*u);
u+=eph->cus*sin2u+eph->cuc*cos2u;
r+=eph->crs*sin2u+eph->crc*cos2u;
i+=eph->cis*sin2u+eph->cic*cos2u;
x=r*cos(u); y=r*sin(u); cosi=cos(i);

基于广播星历给出的六个摄动参数，据此计算出卫星初始三个参数的摄动力改正项，但是这里对于轨道倾角的改正并不彻底，没有调用轨道倾角的变化率参数来换算到此时卫星倾角的实际值。

计算卫星在轨道坐标系中的位置

x=r*cos(u); y=r*sin(u); cosi=cos(i);

在轨道平面直角坐标系中（坐标原点位于地心，X轴指向升交点）卫星的平面直角坐标为：x=rcosu, y=rsinu

计算观测瞬间升交点在协议地球坐标系中的赤经

        O=eph->OMG0+(eph->OMGd-omge)*tk-omge*eph->toes;

利用参考时刻的升交点赤经加上到参考时刻升交点赤经的变化，得到观测时刻的升交点赤经，但是由于地球自转的原因，在地固系的坐标系随着地球旋转，想要得到协议地球坐标系中的升交点赤经，必须减去地球自转所造成的卫星赤经变化，得到观测时刻升交点在协议地球坐标系中的赤经。

计算卫星在协议地球坐标系中的位置

sinO=sin(O); cosO=cos(O);
rs[0]=x*cosO-y*cosi*sinO;
rs[1]=x*sinO+y*cosi*cosO;
rs[2]=y*sin(i);

已知升交点的大地经度以及轨道的平面倾角i后，就可以通过两次旋转方便的求得卫星在地固系中的位置。

计算观测时刻的卫星钟差初值

 tk=timediff(time,eph->toc);
 *dts=eph->f0+eph->f1*tk+eph->f2*tk*tk;

计算观测时刻与参考时刻toc时间差，利用多项式拟合的方式，根据广播星历的af0、af1、af2算出观测时刻的卫星钟差。

计算观测时刻相对论效应改正后的卫星钟差

*dts-=2.0*sqrt(mu*eph->A)*eph->e*sinE/SQR(CLIGHT);

此处考虑到卫星处于高速的运动状态会对卫星钟的准确度产生影响，使用相对论原理，对这项误差进行改正，输出卫星钟差的参数。

计算由于位置和钟不准对定位产生的方差大小

*var=var_uraeph(eph->sva);
static double var_uraeph(int ura)
{
    const double ura_value[]={   
        2.4,3.4,4.85,6.85,9.65,13.65,24.0,48.0,96.0,192.0,384.0,768.0,1536.0,
        3072.0,6144.0
    };
    return ura<0||15<ura?SQR(6144.0):SQR(ura_value[ura]);
}

此处是使用广播星历对应的URA数值，直接估计出方差的大小，如果ura的数值小于0大于15那么直接就是最大值6144的平方作为返回值，如果在0与15之间，那么返回ura_value数组对应值的平方作为输出值。

加入时间的微小改变量，再次计算卫星的位置、钟差

time=timeadd(time,tt);
eph2pos(time,eph,rst,dtst,var);
*svh=eph->svh;

观测时刻加上0.001时间微小变化量，再次计算卫星的位置钟差。

计算卫星的速度、钟漂

 for (i=0;i<3;i++) rs[i+3]=(rst[i]-rs[i])/tt;
 dts[1]=(dtst[0]-dts[0])/tt;

通过相邻时刻下，位置的变化除以时间得到速度，同理利用钟差的变化除以时间的得到钟漂，到此，卫星位置计算全部结束。

参考

使用广播星历的轨道参数进行卫星定轨数学原理可以参考此文章 GPS卫星位置计算