RTKLIB专题学习(四)---单点定位实现初识（一）

RTKLIB专题学习(四)

今天我们来学习RTKLIB中单点定位的调用情况

1.单点定位的核心主程序，我觉得其实是这一部分

该函数完成了每个历元的接收机位置和钟差参数的解算工作，用到的是RTKLIB专题学习(三)—矩阵应用里面的加权最小二乘

/* estimate receiver position ------------------------------------------------*/
static int estpos(const obsd_t* obs, int n, const double* rs, const double* dts,const double* vare, const int* svh, const nav_t* nav,const prcopt_t* opt, sol_t* sol, double* azel, int* vsat,double* resp, char* msg)
{double x[NX] = { 0 },  ave = 0.0, dx[NX], Q[NX * NX], * v, * H, * H1, * var, * ele, *xp,*Pp,* R, * I, * E, * KV, sig, * P, * L, * Norv,deltax=0.0, c = 0.0005, k0 = 0.1, k1 = 2, rr = 0.0;int i, j, k,m, info, stat, nv, ns, num = 0;trace(3, "estpos  : n=%d\n", n);v = mat(n + 4, 1); H = mat(NX, n + 4); H1 = mat(NX, n + 4); var = mat(n + 4, 1); ele = mat(n + 4, 1);  P = eye(NX); Norv = mat(n + 4, 1); R = zeros(n + 4, n + 4); KV = mat(NX, 1);xp = mat(NX, 1);Pp = mat(NX, NX);for (i = 0; i < 3; i++) x[i] = sol->rr[i];for (i = 0; i < MAXITR; i++) {if (sol->num != 1){nv = rescode(i, obs, n, rs, dts, vare, svh, nav, x, opt, v, H, var, azel, vsat, resp,&ns, ele, R);}else{matcpy(Pp,sol->P, NX, NX);matcpy(xp, sol->X, NX, 1);nv = rescode(i, obs, n, rs, dts, vare, svh, nav, xp, opt, v, H, var, azel, vsat, resp,&ns, ele, R);}if (nv < NX) {sprintf(msg, "lack of valid sats ns=%d", nv);break;}if ((norm(dx, NX) < 1E-4)){matcpy(sol->P,Q,  NX, NX);matcpy(sol->X, x, NX, 1);//sol->num += 1;//只有KF和xg-ls时才打开//sol->nm += 1;//只有KF时才打开sol->type = 0;sol->time = timeadd(obs[0].time, -x[3] / CLIGHT);sol->dtr[0] = x[3] / CLIGHT; /* receiver clock bias (s) */sol->dtr[1] = x[4] / CLIGHT; /* GLO-GPS time offset (s) */sol->dtr[2] = x[5] / CLIGHT; /* GAL-GPS time offset (s) */sol->dtr[3] = x[6] / CLIGHT; /* BDS-GPS time offset (s) */sol->dtr[4] = x[7] / CLIGHT; /* IRN-GPS time offset (s) */for (j = 0; j < 6; j++) sol->rr[j] = j < 3 ? x[j] : 0.0;for (j = 0; j < 3; j++) sol->qr[j] = (float)Q[j + j * NX];sol->qr[3] = (float)Q[1];    /* cov xy */sol->qr[4] = (float)Q[2 + NX]; /* cov yz */sol->qr[5] = (float)Q[2];    /* cov zx */sol->ns = (uint8_t)ns;sol->age = sol->ratio = 0.0;/* validate solution */if ((stat = valsol(azel, vsat, n, opt, v, nv, NX, msg))) {sol->stat = opt->sateph == EPHOPT_SBAS ? SOLQ_SBAS : SOLQ_SINGLE;}if (stat) {estvel(obs, n, rs, dts, nav, opt, sol, azel, vsat);}free(v); free(H); free(H1); free(var); free(ele); free(P);  free(R); free(Norv); free(xp); free(Pp);return stat;}if (i>=MAXITR) sprintf(msg,"iteration divergent i=%d",i);free(v); free(H); free(H1); free(var); free(ele); free(P);  free(Norv); free(xp); free(Pp);return 0;
}

2.然后呢，调用estpos的函数为pntpos

/* single-point positioning ----------------------------------------------------
* compute receiver position, velocity, clock bias by single-point positioning
* with pseudorange and doppler observables
* args   : obsd_t *obs      I   observation data
*          int    n         I   number of observation data
*          nav_t  *nav      I   navigation data
*          prcopt_t *opt    I   processing options
*          sol_t  *sol      IO  solution
*          double *azel     IO  azimuth/elevation angle (rad) (NULL: no output)
*          ssat_t *ssat     IO  satellite status              (NULL: no output)
*          char   *msg      O   error message for error exit
* return : status(1:ok,0:error)
*-----------------------------------------------------------------------------*/
extern int pntpos(const obsd_t *obs, int n, const nav_t *nav,const prcopt_t *opt, sol_t *sol, double *azel, ssat_t *ssat,char *msg)
{prcopt_t opt_=*opt;double *rs,*dts,*var,*azel_,*resp,*ele;int i,stat,vsat[MAXOBS]={0},svh[MAXOBS];trace(3,"pntpos  : tobs=%s n=%d\n",time_str(obs[0].time,3),n);sol->stat=SOLQ_NONE;if (n<=0) {strcpy(msg,"no observation data");return 0;}sol->time=obs[0].time;msg[0]='\0';rs = mat(6, n); dts = mat(2, n); var = mat(1, n); azel_ = zeros(2, n); resp = mat(1, n); ele = mat(1,n);if (opt_.mode!=PMODE_SINGLE) { /* for precise positioning */opt_.ionoopt=IONOOPT_BRDC;opt_.tropopt=TROPOPT_SAAS;}/* satellite positons, velocities and clocks */satposs(sol->time,obs,n,nav,opt_.sateph,rs,dts,var,svh);/* estimate receiver position with pseudorange */stat=estpos(obs,n,rs,dts,var,svh,nav,&opt_,sol,azel_,vsat,resp,msg);/* RAIM FDE */if (!stat&&n>=6&&opt->posopt[4]) {stat=raim_fde(obs,n,rs,dts,var,svh,nav,&opt_,sol,azel_,vsat,resp,msg);}/* estimate receiver velocity with Doppler */if (stat) {estvel(obs,n,rs,dts,nav,&opt_,sol,azel_,vsat);}if (azel) {for (i=0;i<n*2;i++) azel[i]=azel_[i];}if (ssat) {for (i=0;i<MAXSAT;i++) {ssat[i].vs=0;ssat[i].azel[0]=ssat[i].azel[1]=0.0;ssat[i].resp[0]=ssat[i].resc[0]=0.0;ssat[i].snr[0]=0;}for (i=0;i<n;i++) {ssat[obs[i].sat-1].azel[0]=azel_[  i*2];ssat[obs[i].sat-1].azel[1]=azel_[1+i*2];ssat[obs[i].sat-1].snr[0]=obs[i].SNR[0];if (!vsat[i]) continue;ssat[obs[i].sat-1].vs=1;ssat[obs[i].sat-1].resp[0]=resp[i];}}free(rs); free(dts); free(var); free(azel_); free(resp);return stat;
}

3.在pntpos中，还调用了satposs卫星位置、速度和钟差函数；raim_fde错误探测和排除函数；以及estvel接收机速度估计函数
4.然后呢，调用pntpos函数的是rtkpos，需要强调的是，rtkpos在精密单点定位中也要调用，在这个函数里面，判定了是要进行单点定位还是精密单点定位

/* precise positioning ---------------------------------------------------------
* input observation data and navigation message, compute rover position by
* precise positioning
* args   : rtk_t *rtk       IO  RTK control/result struct
*            rtk->sol       IO  solution
*                .time      O   solution time
*                .rr[]      IO  rover position/velocity
*                               (I:fixed mode,O:single mode)
*                .dtr[0]    O   receiver clock bias (s)
*                .dtr[1-5]  O   receiver GLO/GAL/BDS/IRN/QZS-GPS time offset (s)
*                .Qr[]      O   rover position covarinace
*                .stat      O   solution status (SOLQ_???)
*                .ns        O   number of valid satellites
*                .age       O   age of differential (s)
*                .ratio     O   ratio factor for ambiguity validation
*            rtk->rb[]      IO  base station position/velocity
*                               (I:relative mode,O:moving-base mode)
*            rtk->nx        I   number of all states
*            rtk->na        I   number of integer states
*            rtk->ns        O   number of valid satellites in use
*            rtk->tt        O   time difference between current and previous (s)
*            rtk->x[]       IO  float states pre-filter and post-filter
*            rtk->P[]       IO  float covariance pre-filter and post-filter
*            rtk->xa[]      O   fixed states after AR
*            rtk->Pa[]      O   fixed covariance after AR
*            rtk->ssat[s]   IO  satellite {s+1} status
*                .sys       O   system (SYS_???)
*                .az   [r]  O   azimuth angle   (rad) (r=0:rover,1:base)
*                .el   [r]  O   elevation angle (rad) (r=0:rover,1:base)
*                .vs   [r]  O   data valid single     (r=0:rover,1:base)
*                .resp [f]  O   freq(f+1) pseudorange residual (m)
*                .resc [f]  O   freq(f+1) carrier-phase residual (m)
*                .vsat [f]  O   freq(f+1) data vaild (0:invalid,1:valid)
*                .fix  [f]  O   freq(f+1) ambiguity flag
*                               (0:nodata,1:float,2:fix,3:hold)
*                .slip [f]  O   freq(f+1) cycle slip flag
*                               (bit8-7:rcv1 LLI, bit6-5:rcv2 LLI,
*                                bit2:parity unknown, bit1:slip)
*                .lock [f]  IO  freq(f+1) carrier lock count
*                .outc [f]  IO  freq(f+1) carrier outage count
*                .slipc[f]  IO  freq(f+1) cycle slip count
*                .rejc [f]  IO  freq(f+1) data reject count
*                .gf        IO  geometry-free phase (L1-L2) (m)
*                .gf2       IO  geometry-free phase (L1-L5) (m)
*            rtk->nfix      IO  number of continuous fixes of ambiguity
*            rtk->neb       IO  bytes of error message buffer
*            rtk->errbuf    IO  error message buffer
*            rtk->tstr      O   time string for debug
*            rtk->opt       I   processing options
*          obsd_t *obs      I   observation data for an epoch
*                               obs[i].rcv=1:rover,2:reference
*                               sorted by receiver and satellte
*          int    n         I   number of observation data
*          nav_t  *nav      I   navigation messages
* return : status (0:no solution,1:valid solution)
* notes  : before calling function, base station position rtk->sol.rb[] should
*          be properly set for relative mode except for moving-baseline
*-----------------------------------------------------------------------------*/
extern int rtkpos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
{prcopt_t *opt=&rtk->opt;sol_t solb={{0}};gtime_t time;int i,nu,nr;char msg[128]="";trace(3,"rtkpos  : time=%s n=%d\n",time_str(obs[0].time,3),n);trace(4,"obs=\n"); traceobs(4,obs,n);/* set base staion position */if (opt->refpos<=POSOPT_RINEX&&opt->mode!=PMODE_SINGLE&&opt->mode!=PMODE_MOVEB) {for (i=0;i<6;i++) rtk->rb[i]=i<3?opt->rb[i]:0.0;}/* count rover/base station observations */for (nu=0;nu   <n&&obs[nu   ].rcv==1;nu++) ;for (nr=0;nu+nr<n&&obs[nu+nr].rcv==2;nr++) ;time=rtk->sol.time; /* previous epoch *//* rover position by single point positioning */if (!pntpos(obs,nu,nav,&rtk->opt,&rtk->sol,NULL,rtk->ssat,msg)) {errmsg(rtk,"point pos error (%s)\n",msg);if (!rtk->opt.dynamics) {outsolstat(rtk);return 0;}}if (time.time != 0) rtk->tt = timediff(rtk->sol.time, time); opt->tt = rtk->tt;/* single point positioning */if (opt->mode==PMODE_SINGLE) {outsolstat(rtk);return 1;}/* suppress output of single solution */if (!opt->outsingle) {rtk->sol.stat=SOLQ_NONE;}/* precise point positioning */if (opt->mode>=PMODE_PPP_KINEMA) {pppos(rtk,obs,nu,nav);outsolstat(rtk);return 1;}/* check number of data of base station and age of differential */if (nr==0) {errmsg(rtk,"no base station observation data for rtk\n");outsolstat(rtk);return 1;}if (opt->mode==PMODE_MOVEB) { /*  moving baseline *//* estimate position/velocity of base station */if (!pntpos(obs+nu,nr,nav,&rtk->opt,&solb,NULL,NULL,msg)) {errmsg(rtk,"base station position error (%s)\n",msg);return 0;}rtk->sol.age=(float)timediff(rtk->sol.time,solb.time);if (fabs(rtk->sol.age)>TTOL_MOVEB) {errmsg(rtk,"time sync error for moving-base (age=%.1f)\n",rtk->sol.age);return 0;}for (i=0;i<6;i++) rtk->rb[i]=solb.rr[i];/* time-synchronized position of base station */for (i=0;i<3;i++) rtk->rb[i]+=rtk->rb[i+3]*rtk->sol.age;}else {rtk->sol.age=(float)timediff(obs[0].time,obs[nu].time);if (fabs(rtk->sol.age)>opt->maxtdiff) {errmsg(rtk,"age of differential error (age=%.1f)\n",rtk->sol.age);outsolstat(rtk);return 1;}}/* relative potitioning */relpos(rtk,obs,nu,nr,nav);outsolstat(rtk);return 1;
}

5.rtkpos中调用了pntpos，即上述单点定位过程中的调用流程；在精密单点定位时，调用pppos；在相对定位时，调用了relpos
6.然后呢，调用rtkpos的是procpos

/* process positioning -------------------------------------------------------*/
static void procpos(FILE *fp, const prcopt_t *popt, const solopt_t *sopt,int mode)
{gtime_t time={0};sol_t sol={{0}};rtk_t rtk;obsd_t obs[MAXOBS*2]; /* for rover and base */double rb[3]={0};int i,nobs,n,solstatic,pri[]={6,1,2,3,4,5,1,6};trace(3,"procpos : mode=%d\n",mode);solstatic=sopt->solstatic&&(popt->mode==PMODE_STATIC||popt->mode==PMODE_PPP_STATIC);rtkinit(&rtk,popt);rtcm_path[0]='\0';while ((nobs=inputobs(obs,rtk.sol.stat,popt))>=0) {/* exclude satellites */for (i=n=0;i<nobs;i++) {if ((satsys(obs[i].sat,NULL)&popt->navsys)&&popt->exsats[obs[i].sat-1]!=1) obs[n++]=obs[i];}if (n<=0) continue;/* carrier-phase bias correction */if (!strstr(popt->pppopt,"-ENA_FCB")) {corr_phase_bias_ssr(obs,n,&navs);}if (!rtkpos(&rtk,obs,n,&navs)) continue;if (mode==0) { /* forward/backward */if (!solstatic) {outsol(fp,&rtk.sol,rtk.rb,sopt);}else if (time.time==0||pri[rtk.sol.stat]<=pri[sol.stat]) {sol=rtk.sol;for (i=0;i<3;i++) rb[i]=rtk.rb[i];if (time.time==0||timediff(rtk.sol.time,time)<0.0) {time=rtk.sol.time;}}}else if (!revs) { /* combined-forward */if (isolf>=nepoch) return;solf[isolf]=rtk.sol;for (i=0;i<3;i++) rbf[i+isolf*3]=rtk.rb[i];isolf++;}else { /* combined-backward */if (isolb>=nepoch) return;solb[isolb]=rtk.sol;for (i=0;i<3;i++) rbb[i+isolb*3]=rtk.rb[i];isolb++;}}if (mode==0&&solstatic&&time.time!=0.0) {sol.time=time;outsol(fp,&sol,rb,sopt);}rtkfree(&rtk);
}

7.在procpos中就是在此处调用rtkpos从而进行参数解算的：if (!rtkpos(&rtk,obs,n,&navs)) continue;
8.然后呢，调用procpos的函数是execses

/* execute processing session ------------------------------------------------*/
static int execses(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt,const solopt_t *sopt, const filopt_t *fopt, int flag,char **infile, const int *index, int n, char *outfile)
{FILE *fp;prcopt_t popt_=*popt;char tracefile[1024],statfile[1024],path[1024],*ext;trace(3,"execses : n=%d outfile=%s\n",n,outfile);/* open debug trace */if (flag&&sopt->trace>0) {if (*outfile) {strcpy(tracefile,outfile);strcat(tracefile,".trace");}else {strcpy(tracefile,fopt->trace);}traceclose();traceopen(tracefile);tracelevel(sopt->trace);}/* read ionosphere data file */if (*fopt->iono&&(ext=strrchr(fopt->iono,'.'))) {if (strlen(ext)==4&&(ext[3]=='i'||ext[3]=='I')) {reppath(fopt->iono,path,ts,"","");readtec(path,&navs,1);}}/* read erp data */if (*fopt->eop) {free(navs.erp.data); navs.erp.data=NULL; navs.erp.n=navs.erp.nmax=0;reppath(fopt->eop,path,ts,"","");if (!readerp(path,&navs.erp)) {showmsg("error : no erp data %s",path);trace(2,"no erp data %s\n",path);}}/* read obs and nav data */if (!readobsnav(ts,te,ti,infile,index,n,&popt_,&obss,&navs,stas)) return 0;/* read dcb parameters */if (*fopt->dcb) {reppath(fopt->dcb,path,ts,"","");readdcb(path,&navs,stas);}/* set antenna paramters */if (popt_.mode!=PMODE_SINGLE) {setpcv(obss.n>0?obss.data[0].time:timeget(),&popt_,&navs,&pcvss,&pcvsr,stas);}/* read ocean tide loading parameters */if (popt_.mode>PMODE_SINGLE&&*fopt->blq) {readotl(&popt_,fopt->blq,stas);}/* rover/reference fixed position */if (popt_.mode==PMODE_FIXED) {if (!antpos(&popt_,1,&obss,&navs,stas,fopt->stapos)) {freeobsnav(&obss,&navs);return 0;}}else if (PMODE_DGPS<=popt_.mode&&popt_.mode<=PMODE_STATIC) {if (!antpos(&popt_,2,&obss,&navs,stas,fopt->stapos)) {freeobsnav(&obss,&navs);return 0;}}/* open solution statistics */if (flag&&sopt->sstat>0) {strcpy(statfile,outfile);strcat(statfile,".stat");rtkclosestat();rtkopenstat(statfile,sopt->sstat);}/* write header to output file */if (flag&&!outhead(outfile,infile,n,&popt_,sopt)) {freeobsnav(&obss,&navs);return 0;}iobsu=iobsr=isbs=revs=aborts=0;if (popt_.mode==PMODE_SINGLE||popt_.soltype==0) {if ((fp=openfile(outfile))) {procpos(fp,&popt_,sopt,0); /* forward */fclose(fp);}}else if (popt_.soltype==1) {if ((fp=openfile(outfile))) {revs=1; iobsu=iobsr=obss.n-1; isbs=sbss.n-1;procpos(fp,&popt_,sopt,0); /* backward */fclose(fp);}}else { /* combined */solf=(sol_t *)malloc(sizeof(sol_t)*nepoch);solb=(sol_t *)malloc(sizeof(sol_t)*nepoch);rbf=(double *)malloc(sizeof(double)*nepoch*3);rbb=(double *)malloc(sizeof(double)*nepoch*3);if (solf&&solb) {isolf=isolb=0;procpos(NULL,&popt_,sopt,1); /* forward */revs=1; iobsu=iobsr=obss.n-1; isbs=sbss.n-1;procpos(NULL,&popt_,sopt,1); /* backward *//* combine forward/backward solutions */if (!aborts&&(fp=openfile(outfile))) {combres(fp,&popt_,sopt);fclose(fp);}}else showmsg("error : memory allocation");free(solf);free(solb);free(rbf);free(rbb);}/* free obs and nav data */freeobsnav(&obss,&navs);return aborts?1:0;
}

篇幅原因，剩下的调用过程会在下一篇进行分析！

RTKLIB专题学习(四)---单点定位实现初识（一）相关推荐

RTKLIB专题学习(四)---单点定位实现初识（二）
RTKLIB专题学习(四)-单点定位实现初识(二) 今天我们来继续学习RTKLIB中单点定位的调用情况,上一篇在这里:RTKLIB专题学习(四)-单点定位实现初识(一) 1.上篇说到了调用procpo ...
RTKLIB专题学习(六)---单点定位应用(二)
RTKLIB专题学习(六)-单点定位应用(二) 上一篇RTKLIB专题学习(六)-单点定位应用(一)我们使用最小二乘以及抗差最小二乘对单点定位进行定位效果分析,发现在存在粗差情况下,抗差最小二乘可以获 ...
RTKLIB专题学习(五)---单点定位实现进阶(一)
RTKLIB专题学习(五) 今天我们一起来了解一下,RTKLIB中的单点定位是如何实现的: 简单来说,就是最小二乘法,但是呢,RTKLIB里面的最小二乘实际上是加权最小二乘,因为他给出了观测值的权(实 ...
RTKLIB专题学习(七)---精密单点定位实现初识（三）
RTKLIB专题学习(七)-精密单点定位实现初识(三) 上两篇我们介绍了RTKLIB中精密单点定位的大致流程,今天我们对照RTKLIB学习手册,来学习相应改正公式和误差源 1.在PPP模式中 RTKL ...
RTKLIB专题学习(十)—电离层改正
RTKLIB专题学习(十)-电离层改正今天我们一起来学习下,RTKLIB中对电离层延迟改正的方法,并了解源码的相应模块,以便对原理有一定的了解 1.电离层延迟改正模型由于存在大量的自由电子和正负离 ...
RTKLIB专题学习(三)---矩阵应用
RTKLIB专题学习(三) 今天我们来进一步学习RTKLIB中矩阵的各种应用 rtkcmn.c : rtklib common functions 1.这是最小二乘法(实际在应用中为加权最小二乘) m ...
IM通讯协议专题学习(八)：金蝶随手记团队的Protobuf应用实践(原理篇)
本文由金蝶随手记技术团队丁同舟分享. 1.引言跟移动端IM中追求数据传输效率.网络流量消耗等需求一样,随手记客户端与服务端交互的过程中,对部分数据的传输大小和效率也有较高的要求,普通的数据格式如 J ...
【响应式编程的思维艺术】（1）Rxjs专题学习计划
[摘要] 请暂时忘掉你的对象,感受一切皆流的世界. 一. 响应式编程响应式编程,也称为流式编程,对于非前端工程师来说,可能并不是一个陌生的名词,它是函数式编程在软件开发中应用的延伸,如果你对函数式编 ...
MySQL学习（一）—初识MySQL
MySQL学习(一)-初识MySQL 害,不知不觉一个半月过去了,时间过得太快了,学校的MySQL已经上到第七周的课程来了,我才想起来要写总结,新学的知识三天不看就忘得差不多了,果然还是要多复习啊,加 ...
Jmeter(总篇): 针对性能测试工具：Jmeter的专题学习
今天因为要提供给其他同事学习的资料,进行整理Jmeter学习资料. 根据我的博客,整理的针对性能测试工具:Jmeter的专题学习,请参考. 第一阶段:组件介绍序号标题链接 1 1.0Jmeter ...

RTKLIB专题学习(四)---单点定位实现初识（一）

RTKLIB专题学习(四)

RTKLIB专题学习(四)---单点定位实现初识（一）相关推荐

最新文章

热门文章