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

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;
}

RTKLIB专题学习(八)—卫星星历和钟差相关推荐

RTKLIB专题学习(十)—电离层改正
RTKLIB专题学习(十)-电离层改正今天我们一起来学习下,RTKLIB中对电离层延迟改正的方法,并了解源码的相应模块,以便对原理有一定的了解 1.电离层延迟改正模型由于存在大量的自由电子和正负离 ...
RTKLIB专题学习(五)---单点定位实现进阶(一)
RTKLIB专题学习(五) 今天我们一起来了解一下,RTKLIB中的单点定位是如何实现的: 简单来说,就是最小二乘法,但是呢,RTKLIB里面的最小二乘实际上是加权最小二乘,因为他给出了观测值的权(实 ...
RTKLIB专题学习(七)---精密单点定位实现初识（三）
RTKLIB专题学习(七)-精密单点定位实现初识(三) 上两篇我们介绍了RTKLIB中精密单点定位的大致流程,今天我们对照RTKLIB学习手册,来学习相应改正公式和误差源 1.在PPP模式中 RTKL ...
RTKLIB专题学习(三)---矩阵应用
RTKLIB专题学习(三) 今天我们来进一步学习RTKLIB中矩阵的各种应用 rtkcmn.c : rtklib common functions 1.这是最小二乘法(实际在应用中为加权最小二乘) m ...
RTKLIB专题学习(四)---单点定位实现初识（二）
RTKLIB专题学习(四)-单点定位实现初识(二) 今天我们来继续学习RTKLIB中单点定位的调用情况,上一篇在这里:RTKLIB专题学习(四)-单点定位实现初识(一) 1.上篇说到了调用procpo ...
RTKLIB专题学习(六)---单点定位应用(二)
RTKLIB专题学习(六)-单点定位应用(二) 上一篇RTKLIB专题学习(六)-单点定位应用(一)我们使用最小二乘以及抗差最小二乘对单点定位进行定位效果分析,发现在存在粗差情况下,抗差最小二乘可以获 ...
IM通讯协议专题学习(八)：金蝶随手记团队的Protobuf应用实践(原理篇)
本文由金蝶随手记技术团队丁同舟分享. 1.引言跟移动端IM中追求数据传输效率.网络流量消耗等需求一样,随手记客户端与服务端交互的过程中,对部分数据的传输大小和效率也有较高的要求,普通的数据格式如 J ...
【python】读取卫星星历（RENIX 3.04）进行卫星位置的计算（北斗卫星专题）
最近的卫星导航数据处理,老师让我们进行卫星位置的计算,从而使用绘图工具进行对卫星星下点的轨迹进行绘图,这里首先的步骤是读取卫星星历数据,计算卫星位置. 这次的课程目标主要是针对北斗卫星,进行对卫星位置 ...
莫比乌斯反演专题学习笔记
莫比乌斯反演专题学习笔记本文记录一些和莫反有关的内容的笔记可能存在诸多谬误,阅读时请谨慎分析若发现文中有谬误,如您愿意,恳请您向我指出,不胜感激! 为什么要学莫比乌斯反演? 解决一类与狄利克雷卷 ...
Python模块EasyGui专题学习
Python模块EasyGui专题学习 1.msgbox(msg,title,ok_button="OK",image="",root=None) 代码 imp ...

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

RTKLIB专题学习(八)—卫星星历和钟差相关推荐

最新文章

热门文章