RTKLIB专题学习—卫星星历和钟差[通俗易懂]-卫星星历分为哪两种

RTKLIB专题学习(八)—卫星星历和钟差[通俗易懂]RTKLIB专题学习(八)—卫星星历和钟差今天我们来学习一下，RTKLIB中使用的卫星星历和钟差情况1.RTKLIB支持GPS、GLONASS、Galileo、QZSS、北斗和SBAS的广播星历和时钟，对于后处理也支持精密星历和时钟(1)GPS、伽利略和QZSS的广播星历和时钟如上公式可知，使用广播星历中参数解算的的是基于卫星天线相位中心的位置，具体公式如下：详细结算流程如下：2.对于GLONASS来说，广播星历直接给出了其位置信息；北斗广播星历也是给出了相应参数，可通_linuxrtklib解析星历

RTKLIB专题学习(八)—卫星星历和钟差

今天我们来学习一下，RTKLIB中使用的卫星星历和钟差情况

1.RTKLIB支持GPS、GLONASS、Galileo、QZSS、北斗和SBAS的广播星历和时钟，对于后处理也支持精密星历和时钟
(1)GPS、伽利略和QZSS的广播星历和时钟

如上公式可知，使用广播星历中参数解算的的是基于卫星天线相位中心的位置，具体公式如下：


详细结算流程如下：








2.对于GLONASS来说，广播星历直接给出了其位置信息；北斗广播星历也是给出了相应参数，可通过与(1)相似的式子算出位置和钟差信息
3.广播星历获取卫星位置和钟差的RTKLIB源码如下：

/* satellite position and clock by broadcast ephemeris -----------------------*/static int ephpos(gtime_t time, gtime_t teph, int sat, const nav_t *nav,int iode, double *rs, double *dts, double *var, int *svh){ eph_t  *eph;geph_t *geph;seph_t *seph;double rst[3],dtst[1],tt=1E-3;int i,sys;trace(4,"ephpos : time=%s sat=%2d iode=%d\n",time_str(time,3),sat,iode);sys=satsys(sat,NULL);*svh=-1;if (sys==SYS_GPS||sys==SYS_GAL||sys==SYS_QZS||sys==SYS_CMP||sys==SYS_IRN) { if (!(eph=seleph(teph,sat,iode,nav))) return 0;eph2pos(time,eph,rs,dts,var);time=timeadd(time,tt);eph2pos(time,eph,rst,dtst,var);*svh=eph->svh;}else if (sys==SYS_GLO) { if (!(geph=selgeph(teph,sat,iode,nav))) return 0;geph2pos(time,geph,rs,dts,var);time=timeadd(time,tt);geph2pos(time,geph,rst,dtst,var);*svh=geph->svh;}else if (sys==SYS_SBS) { if (!(seph=selseph(teph,sat,nav))) return 0;seph2pos(time,seph,rs,dts,var);time=timeadd(time,tt);seph2pos(time,seph,rst,dtst,var);*svh=seph->svh;}else return 0;/* satellite velocity and clock drift by differential approx */for (i=0;i<3;i++) rs[i+3]=(rst[i]-rs[i])/tt;dts[1]=(dtst[0]-dts[0])/tt;return 1;}

希望我今天分享的这篇文章可以帮到您。

其中核心函数为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);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;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));}if (n>=MAX_ITER_KEPLER) { trace(2,"eph2pos: kepler iteration overflow sat=%2d\n",eph->sat);return;}sinE=sin(E); cosE=cos(E);trace(4,"kepler: sat=%2d e=%8.5f n=%2d del=%10.3e\n",eph->sat,eph->e,n,E-Ek);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);/* beidou geo satellite */if (sys==SYS_CMP&&(prn<=5||prn>=59)) {  /* ref [9] table 4-1 */O=eph->OMG0+eph->OMGd*tk-omge*eph->toes;sinO=sin(O); cosO=cos(O);xg=x*cosO-y*cosi*sinO;yg=x*sinO+y*cosi*cosO;zg=y*sin(i);sino=sin(omge*tk); coso=cos(omge*tk);rs[0]= xg*coso+yg*sino*COS_5+zg*sino*SIN_5;rs[1]=-xg*sino+yg*coso*COS_5+zg*coso*SIN_5;rs[2]=-yg*SIN_5+zg*COS_5;}else { 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);}tk=timediff(time,eph->toc);*dts=eph->f0+eph->f1*tk+eph->f2*tk*tk;/* relativity correction */*dts-=2.0*sqrt(mu*eph->A)*eph->e*sinE/SQR(CLIGHT);/* position and clock error variance */*var=var_uraeph(sys,eph->sva);}

可以有上述代码获知，最终计算得到的卫星钟差经过了相对论效应改正
4.精密星历获取卫星位置和钟差的RTKLIB源码如下：

/* satellite position/clock by precise ephemeris/clock ------------------------- * compute satellite position/clock with precise ephemeris/clock * args : gtime_t time I time (gpst) * int sat I satellite number * nav_t *nav I navigation data * int opt I sat postion option * (0: center of mass, 1: antenna phase center) * double *rs O sat position and velocity (ecef) * {x,y,z,vx,vy,vz} (m|m/s) * double *dts O sat clock {bias,drift} (s|s/s) * double *var IO sat position and clock error variance (m) * (NULL: no output) * return : status (1:ok,0:error or data outage) * notes : clock includes relativistic correction but does not contain code bias * before calling the function, nav->peph, nav->ne, nav->pclk and * nav->nc must be set by calling readsp3(), readrnx() or readrnxt() * if precise clocks are not set, clocks in sp3 are used instead *-----------------------------------------------------------------------------*/extern int peph2pos(gtime_t time, int sat, const nav_t *nav, int opt,double *rs, double *dts, double *var){ gtime_t time_tt;double rss[3],rst[3],dtss[1],dtst[1],dant[3]={ 0},vare=0.0,varc=0.0,tt=1E-3;int i;trace(4,"peph2pos: time=%s sat=%2d opt=%d\n",time_str(time,3),sat,opt);if (sat<=0||MAXSAT<sat) return 0;/* satellite position and clock bias */if (!pephpos(time,sat,nav,rss,dtss,&vare,&varc)||!pephclk(time,sat,nav,dtss,&varc)) return 0;time_tt=timeadd(time,tt);if (!pephpos(time_tt,sat,nav,rst,dtst,NULL,NULL)||!pephclk(time_tt,sat,nav,dtst,NULL)) return 0;/* satellite antenna offset correction */if (opt) { satantoff(time,rss,sat,nav,dant);}for (i=0;i<3;i++) { rs[i  ]=rss[i]+dant[i];rs[i+3]=(rst[i]-rss[i])/tt;}/* relativistic effect correction */if (dtss[0]!=0.0) { dts[0]=dtss[0]-2.0*dot(rs,rs+3,3)/CLIGHT/CLIGHT;dts[1]=(dtst[0]-dtss[0])/tt;}else {  /* no precise clock */dts[0]=dts[1]=0.0;}if (var) *var=vare+varc;return 1;}

其中两个主要核心函数为pephpos和pephclk：

/* satellite position by precise ephemeris -----------------------------------*/static int pephpos(gtime_t time, int sat, const nav_t *nav, double *rs,double *dts, double *vare, double *varc){ double t[NMAX+1],p[3][NMAX+1],c[2],*pos,std=0.0,s[3],sinl,cosl;int i,j,k,index;trace(4,"pephpos : time=%s sat=%2d\n",time_str(time,3),sat);rs[0]=rs[1]=rs[2]=dts[0]=0.0;if (nav->ne<NMAX+1||timediff(time,nav->peph[0].time)<-MAXDTE||timediff(time,nav->peph[nav->ne-1].time)>MAXDTE) { trace(3,"no prec ephem %s sat=%2d\n",time_str(time,0),sat);return 0;}/* binary search */for (i=0,j=nav->ne-1;i<j;) { k=(i+j)/2;if (timediff(nav->peph[k].time,time)<0.0) i=k+1; else j=k;}index=i<=0?0:i-1;/* polynomial interpolation for orbit */i=index-(NMAX+1)/2;if (i<0) i=0; else if (i+NMAX>=nav->ne) i=nav->ne-NMAX-1;for (j=0;j<=NMAX;j++) { t[j]=timediff(nav->peph[i+j].time,time);if (norm(nav->peph[i+j].pos[sat-1],3)<=0.0) { trace(3,"prec ephem outage %s sat=%2d\n",time_str(time,0),sat);return 0;}}for (j=0;j<=NMAX;j++) { pos=nav->peph[i+j].pos[sat-1];/* correciton for earh rotation ver.2.4.0 */sinl=sin(OMGE*t[j]);cosl=cos(OMGE*t[j]);p[0][j]=cosl*pos[0]-sinl*pos[1];p[1][j]=sinl*pos[0]+cosl*pos[1];p[2][j]=pos[2];}for (i=0;i<3;i++) { rs[i]=interppol(t,p[i],NMAX+1);}if (vare) { for (i=0;i<3;i++) s[i]=nav->peph[index].std[sat-1][i];std=norm(s,3);/* extrapolation error for orbit */if      (t[0   ]>0.0) std+=EXTERR_EPH*SQR(t[0   ])/2.0;else if (t[NMAX]<0.0) std+=EXTERR_EPH*SQR(t[NMAX])/2.0;*vare=SQR(std);}/* linear interpolation for clock */t[0]=timediff(time,nav->peph[index  ].time);t[1]=timediff(time,nav->peph[index+1].time);c[0]=nav->peph[index  ].pos[sat-1][3];c[1]=nav->peph[index+1].pos[sat-1][3];if (t[0]<=0.0) { if ((dts[0]=c[0])!=0.0) { std=nav->peph[index].std[sat-1][3]*CLIGHT-EXTERR_CLK*t[0];}}else if (t[1]>=0.0) { if ((dts[0]=c[1])!=0.0) { std=nav->peph[index+1].std[sat-1][3]*CLIGHT+EXTERR_CLK*t[1];}}else if (c[0]!=0.0&&c[1]!=0.0) { dts[0]=(c[1]*t[0]-c[0]*t[1])/(t[0]-t[1]);i=t[0]<-t[1]?0:1;std=nav->peph[index+i].std[sat-1][3]+EXTERR_CLK*fabs(t[i]);}else { dts[0]=0.0;}if (varc) *varc=SQR(std);return 1;}

/* satellite clock by precise clock ------------------------------------------*/static int pephclk(gtime_t time, int sat, const nav_t *nav, double *dts,double *varc){ double t[2],c[2],std;int i,j,k,index;trace(4,"pephclk : time=%s sat=%2d\n",time_str(time,3),sat);if (nav->nc<2||timediff(time,nav->pclk[0].time)<-MAXDTE||timediff(time,nav->pclk[nav->nc-1].time)>MAXDTE) { trace(3,"no prec clock %s sat=%2d\n",time_str(time,0),sat);return 1;}/* binary search */for (i=0,j=nav->nc-1;i<j;) { k=(i+j)/2;if (timediff(nav->pclk[k].time,time)<0.0) i=k+1; else j=k;}index=i<=0?0:i-1;/* linear interpolation for clock */t[0]=timediff(time,nav->pclk[index  ].time);t[1]=timediff(time,nav->pclk[index+1].time);c[0]=nav->pclk[index  ].clk[sat-1][0];c[1]=nav->pclk[index+1].clk[sat-1][0];if (t[0]<=0.0) { if ((dts[0]=c[0])==0.0) return 0;std=nav->pclk[index].std[sat-1][0]*CLIGHT-EXTERR_CLK*t[0];}else if (t[1]>=0.0) { if ((dts[0]=c[1])==0.0) return 0;std=nav->pclk[index+1].std[sat-1][0]*CLIGHT+EXTERR_CLK*t[1];}else if (c[0]!=0.0&&c[1]!=0.0) { dts[0]=(c[1]*t[0]-c[0]*t[1])/(t[0]-t[1]);i=t[0]<-t[1]?0:1;std=nav->pclk[index+i].std[sat-1][0]*CLIGHT+EXTERR_CLK*fabs(t[i]);}else { trace(3,"prec clock outage %s sat=%2d\n",time_str(time,0),sat);return 0;}if (varc) *varc=SQR(std);return 1;}