前言：先简单分析main函数到pppos函数流程，最主要目的实现对扩展卡尔曼滤波在PPP中的代码具体分析，对PPP中各项改正不做重要阐述。
在进行PPP代码学习之前，对重要算法的学习至关重要，
可以前往：RTKLIB学习(二)--1、PPP方程和扩展卡尔曼滤波等算法详解-优快云博客
PPP定位主要流程图：
一、流程调用
1、main函数
在main函数这除了输入ppp所需文件，关注点就是关于ppp的处理选项设置了
具体可前往：RTKLIB学习开篇--将命令行程序转为代码调试-优快云博客
ret=postpos(ts,te,tint,0.0,&prcopt,&solopt,&filopt,infile,n,outfile,"","");在main函数中一一掌握上式参数就可以前往postpos函数了。
2、postpos函数
/* open processing session 读取各参数文件*/
if (!openses(popt,sopt,fopt,&navs,&pcvss,&pcvsr)) return -1;读取传入进来的参数设置及参数文件（用于各项改正的数据）传给相应的指针，在rtklib.h中自行找到对应的结构体并学习，后续会一直用到。
if (ts.time!=0&&te.time!=0&&tu>=0.0) {//处理所选历元ts.time是起始历元时间，te.time是结束历元时间，tu默认为0后续解读这是什么
当ts和te不为0时，表明用户手动设置了历元处理起始时间，因此需要进行复杂的时间处理过程。
当起始历元不为0时或起始历元为0时分别进入不同的流程。为了能够多一点了解代码，就从最复杂的开始分析吧！！！
if (ts.time!=0&&te.time!=0&&tu>=0.0) {//处理所选历元
if (timediff(te,ts)<0.0) {//确保时间段正确
showmsg("error : no period");
closeses(&navs,&pcvss,&pcvsr);
return 0;
}timediff()实现求解两个gtime_t数据的时间差,当时间差小于0时终止，在rtklib.h中找到对应定义
typedef struct { /* time struct */
time_t time; /* time (s) expressed by standard time_t */
double sec; /* fraction of second under 1 s */
} gtime_t;time_t实际上为长整型，存储了1970年1月1日0时0分0秒到历元的秒数，sec自然为不足一秒的小数部分。注意gtime_t 不是gps的时间！！！
for (i=0;i<MAXINFILE;i++) {
if (!(ifile[i]=(char *)malloc(1024))) {//分配内存
for (;i>=0;i--) free(ifile[i]);//分配内存失败时，释放内存
closeses(&navs,&pcvss,&pcvsr);
return -1;
}
} if (tu==0.0||tu>86400.0*MAXPRCDAYS) tu=86400.0*MAXPRCDAYS;
settspan(ts,te);
86400对应一天的秒数，据此推断tu应该是起始历元到终止历元的秒数，就叫处理时间段吧。
extern void settspan(gtime_t ts, gtime_t te) {}settspan()在rtkcmn.c文件没有实现任何功能
/* time to gps time ------------------------------------------------------------
* convert gtime_t struct to week and tow in gps time
* args : gtime_t t I gtime_t struct
* int *week IO week number in gps time (NULL: no output)
* return : time of week in gps time (s)
*-----------------------------------------------------------------------------*/
extern double time2gpst(gtime_t t, int *week)
{
gtime_t t0=epoch2time(gpst0);
time_t sec=t.time-t0.time;//gps起始周至起始历元的时间段单位秒
int w=(int)(sec/(86400*7));//所在GPS整周数
if (week) *week=w;
return (double)(sec-w*86400*7)+t.sec;//周内秒+不足1秒的小数部分
} tunit=tu<86400.0?tu:86400.0;//处理历元上限为1天的秒
tss=tunit*(int)floor(time2gpst(ts,&week)/tunit);//周内起始历元时刻已过完整天的秒tss就叫周内整天秒吧
for (i=0;;i++) { /* for each periods */
tts=gpst2time(week,tss+i*tu);
tte=timeadd(tts,tu-DTTOL);当处间段很长（超过一天时），分段设置起始历元和结束历元进行定位，对于一天的定位没影响。
if (timediff(tts,te)>0.0) break;//tts-te
if (timediff(tts,ts)<0.0) tts=ts;
if (timediff(tte,te)>0.0) tte=te;检查时间并校正处理时刻。
for (j=k=nf=0;j<n;j++) {
ext=strrchr(infile[j],'.');//索引到.位置，获取文件后缀
if (ext&&(!strcmp(ext,".rtcm3")||!strcmp(ext,".RTCM3"))) {//实时数据流
strcpy(ifile[nf++],infile[j]);
}
else {
/* include next day precise ephemeris or rinex brdc nav */
ttte=tte;//有精密星历时，增加历元增加一小时的历元
if (ext&&(!strcmp(ext,".sp3")||!strcmp(ext,".SP3")||
!strcmp(ext,".eph")||!strcmp(ext,".EPH"))) {
ttte=timeadd(ttte,3600.0);
}
else if (strstr(infile[j],"brdc")) {//有广播星历时，增加两小时的历元
ttte=timeadd(ttte,7200.0);
}
nf+=reppaths(infile[j],ifile+nf,MAXINFILE-nf,tts,ttte,"","");
} //rtkcmn.c/3070
while (k<nf) index[k++]=j;
if (nf>=MAXINFILE) {
trace(2,"too many input files. trancated\n");
break;
}
}
if (!reppath(outfile,ofile,tts,"","")&&i>0) flag=0;
// replace keywords in file path with date, time, rover and base station id替换文件名内的字符串
/* execute processing session 处理基站，也是进入定位入口*/
stat=execses_b(tts,tte,ti,popt,sopt,fopt,flag,ifile,index,nf,ofile,
rov,base);
if (stat==1) break;
}
for (i=0;i<MAXINFILE;i++) free(ifile[i]);之后进入基站处理函数execse_b
3、execse_b函数
基站处理是针对RTK的，简单过一遍即可。
/* read prec ephemeris and sbas data */
readpreceph(infile,n,popt,&navs,&sbss,&lexs);
for (i=0;i<n;i++) if (strstr(infile[i],"%b")) break;//检索基站ID if (i<n) { /* include base station keywords */当有基站时，进入第一个条件判断，由于此文章关注PPP，不在详细讲解，后续如果有学习RTK的想法时，会在跟进。
else {
stat=execses_r(ts,te,ti,popt,sopt,fopt,flag,infile,index,n,outfile,rov);
}PPP定位则直接进入流动站处理函数execses_r。
4、execses_r函数
/* execute processing session for each rover 分割字符串，识别每个流动站--------------*/
static int execses_r(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,
const char *rov)
{for (i=0;i<n;i++) if (strstr(infile[i],"%r")) break;//检索流动站ID
if (i<n) { /* include rover keywords */同理当有流动站标识符时，进入第一个条件判断，也是RTK定位处理流动站的流程，不详细讲述。
else {
/* execute processing session */
stat=execses(ts,te,ti,popt,sopt,fopt,flag,infile,index,n,outfile);
}PPP直接进入execses函数。
5、execses函数
/* 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);
}通过flag和sopt->trace判读是否打开调试设置trace有6个等级，输出的tracefile会包含代码各种调试信息。
/* read obs and nav data */
if (!readobsnav(ts,te,ti,infile,index,n,&popt_,&obss,&navs,stas)) return 0;
/* 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;
}
}这些都是根据定位模式读取相应的参数文件。说到定位模式如下：
#define PMODE_SINGLE 0 /* positioning mode: single */
#define PMODE_DGPS 1 /* positioning mode: DGPS/DGNSS */
#define PMODE_KINEMA 2 /* positioning mode: kinematic */
#define PMODE_STATIC 3 /* positioning mode: static */
#define PMODE_MOVEB 4 /* positioning mode: moving-base */
#define PMODE_FIXED 5 /* positioning mode: fixed */
#define PMODE_PPP_KINEMA 6 /* positioning mode: PPP-kinemaric */
#define PMODE_PPP_STATIC 7 /* positioning mode: PPP-static */
#define PMODE_PPP_FIXED 8 /* positioning mode: PPP-fixed */如果用到IGS数据PPP一般用到7，如果追求高精度即固定解8也可以。
/* 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;
}根据输出设置输出相应的文件头。后续就是根据popt_.mode和popt_.soltype即定位模式和滤波模式（前向滤波：从头读数据，后向滤波从尾读数据(只适用于后处理)，双向滤波：先进行一次前向滤波再进行一次后向滤波，最后再将两组解结合。结合方式是：每个历元两组进行比较，两组解中若存在固定解，则输出固定解，否则使用smoother函数处理）
参考：GNSS相对定位：RTKLIB之RTK算法流程解读 - 知乎 (zhihu.com)
之后进入procpos函数。
6、procpos函数
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;
//n为一个历元有效观测卫星数if (!rtkpos(&rtk,obs,n,&navs)) continue;正式进入定位模块，后面的条件判断则是一个历元定位结束后根据滤波方向决定下一个读取的历元。
7、rtkpos函数
/* set base staion position */
if (opt->refpos<=3&&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++) ;这些是RTK定位中要处理的流程，这里跳过。nu统计了观测值（卫星）数。
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;
}
} /* precise point positioning */
if (opt->mode>=PMODE_PPP_KINEMA) {
pppos(rtk,obs,nu,nav);
pppoutsolstat(rtk,statlevel,fp_stat);
return 1;
}ppp定位则调用了pppos函数。正式进入PPP定位内容。
二、pppos函数
1、udstate_ppp
/* temporal update of states */
udstate_ppp(rtk,obs,n,nav);首先要知道状态参数：xyz位置，接收机钟差，对流层延迟， phase biases(相位偏差)周跳？这是后续扩展卡尔曼滤波解算的关键部分的。
（1）udpos_ppp
先进行位置参数的更新
/* fixed mode固定解模式 */
if (rtk->opt.mode==PMODE_PPP_FIXED) {
for (i=0;i<3;i++) initx(rtk,rtk->opt.ru[i],1E-8,i);
return;
}
/* initialize position for first epoch */
if (norm(rtk->x,3)<=0.0) {
for (i=0;i<3;i++) initx(rtk,rtk->sol.rr[i],VAR_POS,i);
}
/* static ppp mode 静态PPP模式 */
if (rtk->opt.mode==PMODE_PPP_STATIC) return;
/* kinmatic mode without dynamics */
for (i=0;i<3;i++) {
initx(rtk,rtk->sol.rr[i],VAR_POS,i);
}根据不同的定位模式进行初始化函数initx，只更新前三个位置参数。
int j;
rtk->x[i]=xi;
for (j=0;j<rtk->nx;j++) {
rtk->P[i+j*rtk->nx]=rtk->P[j+i*rtk->nx]=i==j?var:0.0;
}可以看到使用了之前spp定位结果对状态参数及其协方差矩阵进行了更新。
(2)udclk_ppp
原理同位置参数更新
(3)udtrop_ppp
原理同位置参数更新
(4)udbias_ppp
原理同位置参数更新
参数更新详细内容可以参考
RTKLIB专题学习(七)---精密单点定位实现初识（一）_eclipsing satellite block iia-优快云博客
2、satposs
卫星轨道计算部分会大概讲述一下。
/* satellite positions and clocks */
satposs(obs[0].time,obs,n,nav,rtk->opt.sateph,rs,dts,var,svh);计算卫星位置。
在satposs函数中
for (i=0;i<n&&i<MAXOBS;i++) {
for (j=0;j<6;j++) rs [j+i*6]=0.0;
for (j=0;j<2;j++) dts[j+i*2]=0.0;
var[i]=0.0; svh[i]=0;先将卫星位置速度及卫星钟置零。
/* transmission time by satellite clock计算传播时间 */
time[i]=timeadd(obs[i].time,-pr/CLIGHT);/* satellite clock bias by broadcast ephemeris */
if (!ephclk(time[i],teph,obs[i].sat,nav,&dt)) {
trace(2,"no broadcast clock %s sat=%2d\n",time_str(time[i],3),obs[i].sat);
continue;
}（1）ephclk
在此函数中判断不同的卫星系统，并选择星历和计算卫星钟差。感兴趣的自行找到相关函数并解读。
（2）satpos
在此函数中通过精密星历函数peph2pos计算卫星位置。
【1】pephpos
在此函数中进行精密卫星位置的计算，主要用到了轨道多项式插值函数interppol。插值算法用到了内维尔多项式插值，想法是用到一组已知历元的卫星位置通过插值算法，拟合出近似的卫星轨道。属于数值分析的应用，数值分析貌似是研究生的课程，我才大三啊，等我上研究生在细究把。
参考链接：Leon's Blog - 逐次线性插值 (leonfocus.github.io)
【2】pephclk
计算精密卫星钟
【3】satantoff
计算卫星天线相位中心偏移
3、testeclipe
此函数实现排除有遮挡的卫星涉及到计算地球，太阳和月亮相关角的计算，后续会出文章详细讲解算法与代码。此处略过。
4、tidedisp
此处计算地球潮汐引起的位移（固体，海洋，极移），涉及的相关计算暂未学习。后续会出文章详细讲解算法与代码。此处略过。
接着回到pppos函数
xp=mat(rtk->nx,1); Pp=zeros(rtk->nx,rtk->nx);//参数个数和参数协方差
matcpy(xp,rtk->x,rtk->nx,1);//将伪距单点定位状态解赋给xp矩阵
nv=n*rtk->opt.nf*2; v=mat(nv,1); H=mat(rtk->nx,nv); R=mat(nv,nv);nv是观测值个数=有效观测卫星数*载波频率数*2
v是观测向量的残差向量
H是设计矩阵或雅克比矩阵
R是观测值方差阵
for (i=0;i<rtk->opt.niter;i++) {//迭代进入到res_ppp函数完成扩展卡尔曼滤波的预测部分
5、res_ppp
/* earth tides correction */
if (opt->tidecorr) {
tideopt=opt->tidecorr==1?1:7; /* 1:solid, 2:solid+otl+pole */
tidedisp(gpst2utc(obs[0].time),rr,tideopt,&nav->erp,opt->odisp[0],
disp);
for (i=0;i<3;i++) rr[i]+=disp[i];
}再次进行地球潮汐校正
ecef2pos(rr,pos);//坐标系转换 for (i=0;i<n&&i<MAXOBS;i++) {//对H矩阵一列一列赋值从此之后到下面这段代码，进行了各种误差改正获得误差改正数r
/* satellite clock and tropospheric delay */
r+=-CLIGHT*dts[i*2]+dtrp;for (j=0;j<2;j++) { /* for phase and code */
//meas[0]为相位改正距离,m[1]为码改正距离for (k=0;k<nx;k++) H[k+nx*nv]=0.0;//初始值赋0
v[nv]=meas[j]-r;//各项误差改正后的观测值for (k=0;k<3;k++) H[k+nx*nv]=-e[k];//先对位置参数对应系数赋值if (j==0) {//相位
v[nv]-=x[IB(obs[i].sat,opt)];//最终的残差向量
H[IB(obs[i].sat,opt)+nx*nv]=1.0;//最终的雅克比矩阵
}
var[nv]=varerr(obs[i].sat,sys,azel[1+i*2],j,opt)+varm[j]+vare[i]+vart;//残差向量方差
if (j==0) rtk->ssat[sat-1].resc[0]=v[nv];//记录相位残差
else rtk->ssat[sat-1].resp[0]=v[nv];//记录码残差for (i=0;i<nv;i++) for (j=0;j<nv;j++) {
R[i+j*nv]=i==j?var[i]:0.0;//对角线为对应的残差向量方差，其余为0
}先介绍这么多，详细的各项改正分析后续会更新。
之后回到pppos函数
/* measurement update */
matcpy(Pp,rtk->P,rtk->nx,rtk->nx);这段代码复制参数协方差阵后正式进入到扩展卡尔曼滤波的更新部分。
6、filter
扩展卡尔曼滤波的更新部分。
在rtkcmn文件中的EKF代码，通过注释先了解传入参数的定义。
/* kalman filter ---------------------------------------------------------------
* kalman filter state update as follows:
*
* K=P*H*(H'*P*H+R)^-1, xp=x+K*v, Pp=(I-K*H')*P
*
* args : double *x I states vector (n x 1)
* double *P I covariance matrix of states (n x n)
* double *H I transpose of design matrix (n x m)
* double *v I innovation (measurement - model) (m x 1)
* double *R I covariance matrix of measurement error (m x m)
* int n,m I number of states and measurements
* double *xp O states vector after update (n x 1)
* double *Pp O covariance matrix of states after update (n x n)
* return : status (0:ok,<0:error)
* notes : matirix stored by column-major order (fortran convention)
* if state x[i]==0.0, not updates state x[i]/P[i+i*n]定位到filter函数。
extern int filter(double *x, double *P, const double *H, const double *v,
const double *R, int n, int m)ix=imat(n,1); for (i=k=0;i<n;i++) if (x[i]!=0.0&&P[i+i*n]>0.0) ix[k++]=i;存储非零参数和正定协方差矩阵的行索引（将0值排除）。
x_=mat(k,1); xp_=mat(k,1); P_=mat(k,k); Pp_=mat(k,k); H_=mat(k,m);初始化矩阵，并将传入的先验矩阵赋给这些矩阵。
info=filter_(x_,P_,H_,v,R,k,m,xp_,Pp_);//进行矩阵运算 matcpy(Q,R,m,m);
matcpy(xp,x,n,1);
matmul("NN",n,m,n,1.0,P,H,0.0,F); /* Q=H'*P*H+R F=PH*/
matmul("TN",m,m,n,1.0,H,F,1.0,Q);
if (!(info=matinv(Q,m))) {
matmul("NN",n,m,m,1.0,F,Q,0.0,K); /* K=P*H*Q^-1 卡尔曼增益矩阵*/
matmul("NN",n,1,m,1.0,K,v,1.0,xp); /* xp=x+K*v */
matmul("NT",n,n,m,-1.0,K,H,1.0,I); /* Pp=(I-K*H')*P 更新误差协方差矩阵*/
matmul("NN",n,n,n,1.0,I,P,0.0,Pp);
}这些就是按照公式进行解算，还有疑问的可以前往：
RTKLIB学习(二)--1、PPP方程和扩展卡尔曼滤波等算法详解-优快云博客
for (i=0;i<k;i++) {
x[ix[i]]=xp_[i];
for (j=0;j<k;j++) P[ix[i]+ix[j]*n]=Pp_[i+j*k];
}filter_函数是在传入的xp_和Pp_矩阵上对原矩阵进行计算。因此将结果赋给后验矩阵参数x,P。至此，扩展卡尔曼滤波的代码与算法分析结束。
7、pppamb
进行模糊度固定流程。
/* resolve integer ambiguity for ppp -----------------------------------------*/
extern int pppamb(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav,
const double *azel)
{
double elmask;
int i,j,m=0,stat=0,*NW,*sat1,*sat2;
if (n<=0||rtk->opt.ionoopt!=IONOOPT_IFLC||rtk->opt.nf<2) return 0;//无电离层结束
trace(3,"pppamb: time=%s n=%d\n",time_str(obs[0].time,0),n);
//elmask用于确定地理高度角el的阈值
elmask=rtk->opt.elmaskar>0.0?rtk->opt.elmaskar:rtk->opt.elmin;
sat1=imat(n*n,1); sat2=imat(n*n,1); NW=imat(n*n,1);
/* average LC */
average_LC(rtk,obs,n,nav,azel);
/* fix wide-lane ambiguity 固定宽巷模糊度*/
for (i=0;i<n-1;i++) for (j=i+1;j<n;j++) {
//判断每个卫星的条件是否符合条件
if (!rtk->ssat[obs[i].sat-1].vsat[0]||
!rtk->ssat[obs[j].sat-1].vsat[0]||
azel[1+i*2]<elmask||azel[1+j*2]<elmask) continue;
#if 0
/* test already fixed */
if (rtk->ambc[obs[i].sat-1].flags[obs[j].sat-1]&&
rtk->ambc[obs[j].sat-1].flags[obs[i].sat-1]) continue;
#endif
sat1[m]=obs[i].sat;//存储有效卫星的序号
sat2[m]=obs[j].sat;
if (fix_amb_WL(rtk,nav,sat1[m],sat2[m],NW+m)) m++;
}
/* fix narrow-lane ambiguity窄巷模糊度 */
if (rtk->opt.modear==ARMODE_PPPAR) {
stat=fix_amb_ROUND(rtk,sat1,sat2,NW,m);//通过取整
}
else if (rtk->opt.modear==ARMODE_PPPAR_ILS) {
stat=fix_amb_ILS(rtk,sat1,sat2,NW,m);//通过lambda固定
}average_LC()
对PPP中的三频载波LC进行均值计算，并考虑了测量噪声方差、模糊度的历元间隔和周跳的处理。这些信息将用于PPP状态估计中对模糊度参数的约束
for (i=0;i<n;i++) {//循环处理每个卫星
sat=obs[i].sat;
//通过观测卫星的仰角 azel，判断是否可见
if (azel[1+2*i]<rtk->opt.elmin) continue;//仰角小于设定的最小可见仰角，则忽略该卫星
if (satsys(sat,NULL)!=SYS_GPS) continue;//只处理 GPS 系统的卫星
//计算其三频载波 LC。这里使用的 L_LC 和 P_LC 函数是对三频观测值进行 LC 计算的函数
/* triple-freq carrier and code LC (m) */
LC1=L_LC(1,-1, 0,obs[i].L)-P_LC(1,1,0,obs[i].P);
LC2=L_LC(0, 1,-1,obs[i].L)-P_LC(0,1,1,obs[i].P);
LC3=L_LC(1,-6, 5,obs[i].L)-P_LC(1,1,0,obs[i].P);
sig=sqrt(SQR(rtk->opt.err[1])+SQR(rtk->opt.err[2]/sin(azel[1+2*i])));
/* measurement noise variance (m) 计算三频LC的测量噪声方差*/
var1=var_LC(1,1,0,sig*rtk->opt.eratio[0]);
var2=var_LC(0,1,1,sig*rtk->opt.eratio[0]);
var3=var_LC(1,1,0,sig*rtk->opt.eratio[0]);
amb=rtk->ambc+sat-1;//获取卫星的模糊度信息
//如果该卫星发生周跳或者观测到的历元与上一个历元的时间间隔超过阈值MIN_ARC_GAP，则重新初始化卫星的模糊度信息。
if (rtk->ssat[sat-1].slip[0]||rtk->ssat[sat-1].slip[1]||
rtk->ssat[sat-1].slip[2]||amb->n[0]==0.0||
fabs(timediff(amb->epoch[0],obs[0].time))>MIN_ARC_GAP) {
amb->n[0]=amb->n[1]=amb->n[2]=0.0;
amb->LC[0]=amb->LC[1]=amb->LC[2]=0.0;
amb->LCv[0]=amb->LCv[1]=amb->LCv[2]=0.0;
amb->fixcnt=0;
for (j=0;j<MAXSAT;j++) amb->flags[j]=0;
}//进行均值计算，根据历元间隔以及是否发生周跳，更新模糊度的历元数量、LC 的均值和 LC 的方差
/* averaging */
if (LC1) {
amb->n[0]+=1.0;
amb->LC [0]+=(LC1 -amb->LC [0])/amb->n[0];
amb->LCv[0]+=(var1-amb->LCv[0])/amb->n[0];
}
if (LC2) {
amb->n[1]+=1.0;
amb->LC [1]+=(LC2 -amb->LC [1])/amb->n[1];
amb->LCv[1]+=(var2-amb->LCv[1])/amb->n[1];
}
if (LC3) {
amb->n[2]+=1.0;
amb->LC [2]+=(LC3 -amb->LC [2])/amb->n[2];
amb->LCv[2]+=(var3-amb->LCv[2])/amb->n[2];
}
amb->epoch[0]=obs[0].time;//更新卫星的历元信息，记录当前历元的时间到 amb->epoch[0]
}fix_amb_WL()
计算宽巷模糊度
static int fix_amb_WL(rtk_t *rtk, const nav_t *nav, int sat1, int sat2, int *NW)
{
ambc_t *amb1,*amb2;
double BW,vW,lam_WL=lam_LC(1,-1,0);
//获取两颗卫星对应的宽巷模糊度信息
amb1=rtk->ambc+sat1-1;
amb2=rtk->ambc+sat2-1;
if (!amb1->n[0]||!amb2->n[0]) return 0;//如果其中一颗卫星的宽巷模糊度信息不存在，则返回0
/* wide-lane ambiguity */
#ifndef REV_WL_FCB //根据两颗卫星的宽巷LC和宽巷波长lam_LC，计算宽巷模糊度
BW=(amb1->LC[0]-amb2->LC[0])/lam_WL+nav->wlbias[sat1-1]-nav->wlbias[sat2-1];
#else
BW=(amb1->LC[0]-amb2->LC[0])/lam_WL-nav->wlbias[sat1-1]+nav->wlbias[sat2-1];
#endif
*NW=ROUND(BW);//四舍五入得到整数值
//计算宽巷模糊度的方差 vW
/* variance of wide-lane ambiguity */
vW=(amb1->LCv[0]/amb1->n[0]+amb2->LCv[0]/amb2->n[0])/SQR(lam_WL);
//判断计算得到的整数宽巷模糊度是否在一定的阈值范围内，同时考虑其方差
/* validation of integer wide-lane ambigyity */
return fabs(*NW-BW)<=rtk->opt.thresar[2]&&
conffunc(*NW,BW,sqrt(vW))>=rtk->opt.thresar[1];//置信度
}fix_amb_ROUND()
计算窄巷模糊度
for (i=0;i<n;i++) {//遍历每个卫星
j=IB(sat1[i],&rtk->opt);
k=IB(sat2[i],&rtk->opt);
/* narrow-lane ambiguity */
B1=(rtk->x[j]-rtk->x[k]+C2*lam2*NW[i])/lam_NL;//估计
N1=ROUND(B1);//取整
//通过状态协方差矩阵rtk->P中的对应元素，计算窄巷模糊度的方差v1
/* variance of narrow-lane ambiguity */
var[m]=rtk->P[j+j*rtk->nx]+rtk->P[k+k*rtk->nx]-2.0*rtk->P[j+k*rtk->nx];
v1=var[m]/SQR(lam_NL);
/* validation of narrow-lane ambiguity */
if (fabs(N1-B1)>rtk->opt.thresar[2]||
conffunc(N1,B1,sqrt(v1))<rtk->opt.thresar[1]) {
continue;
}判断计算得到的整数窄巷模糊度是否在一定的阈值范围内，同时考虑其方差。验证的条件为整数窄巷模糊度的绝对值与计算得到的窄巷模糊度之差小于 rtk->opt.thresar[2]，并且通过 conffunc 函数计算的置信度大于 rtk->opt.thresar[1]。
/* iono-free ambiguity (m) 使用窄巷模糊度计算无电离层模糊度BC*/
BC=C1*lam1*N1+C2*lam2*(N1-NW[i]);
/* check residuals 检查残差v是否大于阈值THRES_RES*/
v=rtk->ssat[sat1[i]-1].resc[0]-rtk->ssat[sat2[i]-1].resc[0];
vc=v+(BC-(rtk->x[j]-rtk->x[k]));//残差是通过观测值残差v加上无电离层模糊度与状态向量差值之和
/* select fixed ambiguities by dependancy check */
m=sel_amb(sat1,sat2,NC,var,m);//通过模糊度的依赖性检查，选择最终固定的模糊度sel_amb（）
/* select fixed ambiguities 从一组卫星对中选择线性独立的模糊度对-------------------------------*/
static int sel_amb(int *sat1, int *sat2, double *N, double *var, int n)
{
int i,j,flgs[MAXSAT]={0},max_flg=0;//flgs用于标记卫星对的线性依赖性，max_flg 用于跟踪当前最大的依赖性标记
/* sort by variance */
for (i=0;i<n;i++) for (j=1;j<n-i;j++) {//对给定的一组卫星对按照其方差从小到大进行排序
if (var[j]>=var[j-1]) continue;
SWAP_I(sat1[j],sat1[j-1]);//排序后，sat1、sat2、N和var数组中的元素将按方差递增的顺序排列
SWAP_I(sat2[j],sat2[j-1]);
SWAP_D(N[j],N[j-1]);
SWAP_D(var[j],var[j-1]);
}//遍历排序后的卫星对，判断每对卫星是否与前面已选卫星对线性独立
/* select linearly independent satellite pair */
for (i=j=0;i<n;i++) {
if (!is_depend(sat1[i],sat2[i],flgs,&max_flg)) continue;//判断两对卫星对是否线性依赖
sat1[j]=sat1[i];//如果不依赖，则将该卫星对的信息保留在输出数组中
sat2[j]=sat2[i];
N[j]=N[i];
var[j++]=var[i];
}
return j;//返回成功选择的线性独立卫星对数目
}fix_sol()
/* constraints to fixed ambiguities */
for (i=0;i<n;i++) {
j=IB(sat1[i],&rtk->opt);
k=IB(sat2[i],&rtk->opt);
v[i]=NC[i]-(rtk->x[j]-rtk->x[k]);
H[j+i*rtk->nx]= 1.0;
H[k+i*rtk->nx]=-1.0;
R[i+i*n]=SQR(CONST_AMB);
}对于每一个固定的模糊度对，构建一个约束方程。这个方程的形式为 NC[i] = x[j] - x[k]，其中 NC[i] 是宽巷或窄巷模糊度的估计值，x[j] 和 x[k] 分别是卫星对应的状态向量元素。这些方程被组合成一个形如 H * x = v 的约束方程系统。
/* update states with constraints对约束方程进行滤波，以更新状态向量x和协方差矩阵 P */
if ((info=filter(rtk->x,rtk->P,H,v,R,rtk->nx,n))) {
trace(1,"filter error (info=%d)\n",info);
free(v); free(H); free(R);
return 0;
}
/* set solution 将更新后的状态向量中的固定模糊度的估计值和协方差矩阵存储到rtk->xa和rtk->Pa 中*/
for (i=0;i<rtk->na;i++) {
rtk->xa[i]=rtk->x[i];
for (j=0;j<rtk->na;j++) {
rtk->Pa[i+j*rtk->na]=rtk->Pa[j+i*rtk->na]=rtk->P[i+j*rtk->nx];
}
}
/* set flags将相应卫星对的标志位设置为1，表示这些卫星对的模糊度已经被成功固定 */
for (i=0;i<n;i++) {
rtk->ambc[sat1[i]-1].flags[sat2[i]-1]=1;
rtk->ambc[sat2[i]-1].flags[sat1[i]-1]=1;
}fix_amb_ILS()
通过整周模糊度转换为无电离层组合模糊度，并使用整周模糊度和非整周模糊度的组合来固定模糊度
for (i=0;i<n;i++) {//遍历每一个卫星
/* check linear independency 检查卫星对是否线性独立*/
if (!is_depend(sat1[i],sat2[i],flgs,&max_flg)) continue;
j=IB(sat1[i],&rtk->opt);
k=IB(sat2[i],&rtk->opt);
/* float narrow-lane ambiguity (cycle) */
B1[m]=(rtk->x[j]-rtk->x[k]+C2*lam2*NW[i])/lam_NL;//计算浮点数窄巷模糊度 B1
N1[m]=ROUND(B1[m]);//四合五入
/* validation of narrow-lane ambiguity比较整数值N1和浮点数值B1的差异，判断是否满足窄巷模糊度的验证条件 */
if (fabs(N1[m]-B1[m])>rtk->opt.thresar[2]) continue;
/* narrow-lane ambiguity transformation matrix构建窄巷模糊度到无电离层组合模糊度的转换矩阵D*/
D[j+m*rtk->nx]= 1.0/lam_NL;
D[k+m*rtk->nx]=-1.0/lam_NL;
sat1[m]=sat1[i];
sat2[m]=sat2[i];
NW[m++]=NW[i];
}
if (m<3) return 0;
/* covariance of narrow-lane ambiguities窄巷模糊度协方差*/
matmul("TN",m,rtk->nx,rtk->nx,1.0,D,rtk->P,0.0,E);
matmul("NN",m,m,rtk->nx,1.0,E,D,0.0,Q);
/* integer least square通过调用lambda函数，使用整数最小二乘法计算整数窄巷模糊度的估计值，并得到相应的协方差矩阵Q*/
if ((info=lambda(m,2,B1,Q,N1,s))) {
trace(2,"lambda error: info=%d\n",info);
return 0;
}
if (s[0]<=0.0) return 0;
rtk->sol.ratio=(float)(MIN(s[1]/s[0],999.9));
//比较整数窄巷模糊度的协方差矩阵的比率，判断是否满足验证条件
/* varidation by ratio-test比率检验方差 */
if (rtk->opt.thresar[0]>0.0&&rtk->sol.ratio<rtk->opt.thresar[0]) {
trace(2,"varidation error: n=%2d ratio=%8.3f\n",m,rtk->sol.ratio);
return 0;
}
trace(2,"varidation ok: %s n=%2d ratio=%8.3f\n",time_str(rtk->sol.time,0),m,
rtk->sol.ratio);
/* narrow-lane to iono-free ambiguity */
for (i=0;i<m;i++) {//将整数窄巷模糊度和非整数窄巷模糊度转换为无电离层组合模糊度
NC[i]=C1*lam1*N1[i]+C2*lam2*(N1[i]-NW[i]);
}
/* fixed solution固定解 */
stat=fix_sol(rtk,sat1,sat2,NC,m);lambda（）
/* LD factorization 对浮点协方差阵进行LD分解*/
if (!(info=LD(n,Q,L,D))) {
/* lambda reduction 降相关*/
reduction(n,L,D,Z);
matmul("TN",n,1,n,1.0,Z,a,0.0,z); /* z=Z'*a */// z变换，将双差模糊度进行变换
/* mlambda search *///调用search(),mlambda search，结果存储在E和s中（整数解）
if (!(info=search(n,m,L,D,z,E,s))) {
/* 将在新空间中固定的模糊度逆变换回双差模糊度空间中 */
info=solve("T",Z,E,n,m,F); /* F=Z'\E */
}//逆Z变换，将在新空间中固定的模糊度逆变换回双差模糊度空间中，存储在F中
}reduction()
int i,j,k;
double del;
//这里的调序变换类似插入排序的思路？
j=n-2; k=n-2;//由于第k+1，k+2，...，n-2列都进行过降相关并且没有被上一次调序变换影响，
while (j>=0) {//因此只需对第0,1，...，k-1，k列进行降相关
if (j<=k) for (i=j+1;i<n;i++) gauss(n,L,Z,i,j);//从最后一列开始，各列非对角线元素从上往下依次降相关
del=D[j]+L[j+1+j*n]*L[j+1+j*n]*D[j+1];
if (del+1E-6<D[j+1]) { /* compared considering numerical error 检验条件，若不满足检验条件则开始进行调序变换*/
perm(n,L,D,j,del,Z);//调序变换
k=j; j=n-2;//完成调序变换后重新从最后一列开始进行降相关及排序，k记录最后一次进行过调序变换的列序号
}
else j--;
}search()
int i,j,k,c,nn=0,imax=0;
double newdist,maxdist=1E99,y;//maxdist，当前超椭圆半径
double *S=zeros(n,n),*dist=mat(n,1),*zb=mat(n,1),*z=mat(n,1),*step=mat(n,1);
k=n-1; dist[k]=0.0; //k表示当前层，从最后一层（n-1）开始计算
zb[k]=zs[k];
z[k]=ROUND(zb[k]); y=zb[k]-z[k]; step[k]=SGN(y);//四舍五入取整；取整后的数与未取整的数作差；step记录z[k]是四舍还是五入
for (c=0;c<LOOPMAX;c++) {
newdist=dist[k]+y*y/D[k];
if (newdist<maxdist) {//如果当前累积目标函数计算值小于当前超椭圆半径
if (k!=0) {//若还未计算至第一层，继续计算累积目标函数值
dist[--k]=newdist;//记录下当前层的累积目标函数值，dist[k]表示了第k,k+1,...,n-1层的目标函数计算和
for (i=0;i<=k;i++)
S[k+i*n]=S[k+1+i*n]+(z[k+1]-zb[k+1])*L[k+1+i*n];//计算Zk，即第k个整数模糊度参数的备选组的中心
zb[k]=zs[k]+S[k+k*n];
z[k]=ROUND(zb[k]); y=zb[k]-z[k]; step[k]=SGN(y);//四舍五入取整；取整后的数与未取整的数作差；记录是四舍还是五入
}
else {//若已经计算至第一层，意味着所有层的累积目标函数值计算完毕
//nn为当前候选解数，m为我们需要的固定解数，这里为2，表示需要一个最优解及一个次优解
//s记录候选解的目标函数值，imax记录之前候选解中的最大目标函数值的坐标
if (nn<m) {//若候选解数还没满
if (nn==0||newdist>s[imax]) imax=nn;//若当前解的目标函数值比之前最大的目标函数值都大，那么更新imax使s[imax]指向当前解中具有的最大目标函数值
for (i=0;i<n;i++) zn[i+nn*n]=z[i];//zn存放所有候选解
s[nn++]=newdist;//s记录当前目标函数值newdist，并加加当前候选解数nn
}
else {//若候选解数已满（即当前zn中已经存了2个候选解）
if (newdist<s[imax]) { //若 当前解的目标函数值 比 s中的最大目标函数值小
for (i=0;i<n;i++) zn[i+imax*n]=z[i];//用当前解替换zn中具有较大目标函数值的解
s[imax]=newdist;//用当前解的目标函数值替换s中的最大目标函数值
for (i=imax=0;i<m;i++) if (s[imax]<s[i]) imax=i;//更新imax保证imax始终指向s中的最大目标函数值
}
maxdist=s[imax];//用当前最大的目标函数值更新超椭圆半径
}//在第一层，取下一个有效的整数模糊度参数进行计算（若zb为5.3，则z取值顺序为5,6,4,7，...）
z[0]+=step[0]; y=zb[0]-z[0]; step[0]=-step[0]-SGN(step[0]);
}
}
else {//如果当前累积目标函数计算值大于当前超椭圆半径
if (k==n-1) break;//如果当前层为第n-1层，意味着后续目标函数各项的计算都会超出超椭圆半径，因此终止搜索
else {//若当前层不是第n-1层
k++;//退后一层，即从第k层退到第k+1层
z[k]+=step[k]; y=zb[k]-z[k]; step[k]=-step[k]-SGN(step[k]);//计算退后一层后，当前层的下一个有效备选解
}
}// 对s中的目标函数值及zn中的候选解进行排序（以s中目标函数值为排序标准，进行升序排序）
}// RTKLIB中最终可以得到一个最优解一个次优解，存在zn中，两解对应的目标函数值，存在s中
for (i=0;i<m-1;i++) { /* sort by s */
for (j=i+1;j<m;j++) {
if (s[i]<s[j]) continue;
SWAP(s[i],s[j]);
for (k=0;k<n;k++) SWAP(zn[k+i*n],zn[k+j*n]);
}
}三、总结
未完，敬请期待后续更新！！
上传一位宝藏博主，找到他的他的Github文档，你会发现意料不到的惊喜
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
