9.2 rtkpos(): 算法入口
本章节，我们主要关注 rtkpos() 中的 relpos() 函数，其中的 satposs() 可以参考第 8 章，udstate() 将在 9.3 节解析，zdres() 和 ddres() 将在 9.4 节进行解析，filter() 函数将在 9.5 节解析，manage_amb_LAMBDA()、resamb_LAMBDA() 和 holdamb() 等与模糊度相关的函数将在 9.6 节解析。
本章主要解析 relpos() 及其调用的部分函数（上述提及的函数除外）。
9.2.1 rtkpos()：单历元 GNSS 定位解算
图9.2-1 rtkpos() 主要组成
1. 参数列表
c
/* args */
rtk_t *rtk IO rtk 控制结构体
const obsd_t *obs I 原始观测数据 obs
int n I 观测数据数目
const nav_t *nav I 导航电文信息
/* return */
int status - 0：无解，1：有效解2. 执行流程
- 设置 RTK 基站坐标
rtk->rb：其中速度设置为 0.0，以后处理为例，基站坐标在execses()函数内调用antpos()函数根据配置选项按照不同的方式获取位置坐标：postype=POSOPT_SINGLE：调用avepos()利用基站的观测文件计算其 SPP 定位结果平均值作为基站的坐标（平均的时间由 ant2-maxaveep 参数进行控制）；postype=POSOPT_FILE：调用getstapos()从 pos 文件中读取基站坐标；postype=POSOPT_RINEX：从 rinex 头文件中获取位置数据。头文件中的测站数据经过读取后已存到stas中；
- 统计 obs 数组：统计基站 obs 个数
nu，流动站 obs 个数nr； - 伪距单点定位：调用
pntpos()计算流动站坐标，作为 Kalman 滤波的坐标初值，如果由于流动站 SPP 定位结果坐标误差过大等原因导致的 SPP 无解，则不进行 RTK 运算，当前历元无解； - 计算前后时间差：计算当前历元和上一历元时间差
rtk->tt（作为时间更新时的 delta_t）； - 单点定位模式：直接输出此前
pntpos()计算得到的坐标； - out-outsingle 选项处理：如果开启了 out-outsingle 选项，那么在 DGNSS/浮点解/固定解/PPP 算法信号中断时（outage），将输出单点定位（SPP）结果。否则将抑制单点解的输出（
rtk->sol.stat=SOLQ_NONE）； - PPP 定位：如果是 PPP 模式，调用
pppos()解算，输出结果； - Moving-Base 模式：
- 调用
pntpos()，传入obs+nu基站观测值计算基站坐标solb； - 计算差分龄期
rtk->sol.age； - 把
solb.rr赋值给rtk->rb； - 时间同步：位置+=对应速度*差分龄期。
- 调用
- 非Moving-Base 模式：计算差分龄期
rtk->sol.age，等于第一个流动站观测值时间obs[0].time减去第一个基站观测值时间obs[nu].time； - 相对定位：调用
relpos()进行 RTK 解算，包含 Kinematic 和 DGNSS 结果； - 解状态输出：调用
outsolstat()输出解算结果，需要开启 out-outstat 选项，它通过控制一个statlevel全局变量来实现。
3. 注意事项
- 在 demo5 b34k 版本代码中，
rtkpos()函数并没有再每个历元都计算一次pntpos()，而是在状态初始化或重置的时候才运行pntpos()，这样可以减少不必要的运算。cif (rtk->P[0]==0||rtk->P[0]>STD_PREC_VAR_THRESH) { if (!pntpos(obs,nu,nav,&rtk->opt,&rtk->sol,NULL,rtk->ssat,msg)) { ... } }
4. 源码注释
点击查看代码
c
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);
// 设置 RTK 基站坐标，以后处理为例，基站坐标在 execses 函数内已经计算了，其中速度设为 0.0
// 这里将配置结构体 opt 中基站的坐标赋值给 rtk 结构体
/* 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; // opt 内基站坐标赋值给 rtk->rb, 速度设为 0.0
}
// 统计基站 obs 个数 nu，流动站 obs 个数 nr，可用于后面判断是否满足差分条件
/* 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 */
// 利用观测值及星历计算流动站的 SPP 定位结果，作为 Kalman 滤波的近似坐标。
// 如果由于流动站 SPP 定位结果坐标误差过大等原因导致的 SPP 无解，则不进行 RTK 运算，当前历元无解。
/* 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;
}
}
// 计算当前历元和上一历元时间差 rtk->tt，rtk->sol.time 是当前历元时间，time 是上一历元时间
if (time.time!=0) rtk->tt=timediff(rtk->sol.time,time);
/* single point positioning */
if (opt->mode==PMODE_SINGLE) { // 单点定位模式直接输出刚刚 SPP 算的坐标
outsolstat(rtk);
return 1;
}
// 如果不是单点模式，抑制单点解的输出（rtk->sol.stat 设为 SOLQ_NONE）
/* 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;
}
// 动基线的基站需要进行实时解算，这也是为什么 rtkpos 初始化时，动基线不参与基站的处理，
// 与其他差分定位相比，动基线要求更严格的实践同步。
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); // 计算差分龄期 rtk->sol.age
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]; // 把 solb.rr 赋值给 rtk->rb
// 时间同步
/* 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;
}9.2.2 relpos()：相对定位
relpos() 是相对定位的主要函数，适用于 DGPS/DGNSS、Kinematic、Static、Moving-Base、Fixed 5 种模式。
图9.2-2 relpos() 函数调用
1. 参数列表
c
/* args */
rtk_t *rtk IO rtk 控制/解算结果结构体
obsd_t *obs I 卫星原始观测数据
int nu I 接收机观测数据数目
int nr I 基站观测数据数目
nav_t *nav I 导航数据
/* return */
int status - 1:ok, 0:error2. 执行流程
- 计算差分龄期：计算流动站、基站间差分龄期
dt：timediff(obs[0].time,obs[nu].time)； - 卫星轨道计算：调用
satposs()计算当前历元下，各卫星的位置速度rs、钟差dts； - 基站非差残差：调用
zdres()计算基站的各卫星观测值的非差残差（观测值-计算值）及卫星的高度角、方位角、卫地距等； - 基站插值处理：后处理中可能会调用
intpres()进行基站残差的插值处理； - 共视卫星选取：调用
selsat()选择基站和流动站之间的共视卫星，进行 RTK 算法时只需要基线间的同步观测卫星，返回共同观测的卫星个数ns，输出卫星号列表sat，在接收机观测值中的 index 值列表iu和在基站观测值中的 index 值列表ir - Kalman 时间更新：调用
udstate()更新状态值rtk->x及其误差协方差rtk->P； - 设置迭代次数
niter，动基线加2次 ； - for 循环迭代量测更新
niter次：- 流动站非差残差：调用
zdres()计算流动站的各位卫星观测值非差残差（观测值-计算值）及卫星的高度角、方位角、卫星矢量等； - 构建双差观测方程：调用
ddres()计算双差残差v，雅可比矩阵H、双差误差协方差矩阵R； - Kalman 量测更新：调用
filter()，计算增益矩阵 K 并进行量测更新，得到 EKF 浮点解。
- 流动站非差残差：调用
- 计算验后残差：量测更新完成，再次调用
zdres()和ddres()计算双差相位/伪距残差，调用valpos()进行浮点解有效性验证；若通过则更新rtk->x以及rtk->P，并更新模糊度控制结构体； - 模糊度的固定：调用
resamb_LAMBDA()，利用 lambda 算法固定模糊度。- 模糊度解算成功，调用
zdres()和ddres()根据固定结果计算残差和协方差，并进行调用valpos()校验； - 固定解验证有效，若为 hold 模式，需要存模糊度信息调用
holdamb()。
- 模糊度解算成功，调用
- 保存中间结果，位置，速度，方差，到
sol.rr、sol.qr、sol.qv，固定解数据在rtk->xa、rtk->Pa，浮点解数据在rtk->x、rtk->P； - 存储当前历元载波信息：
rtk->ssat[sat[i]-1].pt存时间、rtk->ssat[sat[i]-1].ph存载波相位观测值，供下次使用（供detslp_dop使用）； - 存储当前历元 SNR 信息：
rtk->ssat[sat[i]-1].snr； - 存储卫星的模糊度固定信息：固定状态信息
rtk->ssat[i].fix[j]及周跳信息rtk->ssat[i].slipc[j]。
3. 源码注释
点击查看代码
c
static int relpos(rtk_t *rtk, const obsd_t *obs, int nu, int nr,
const nav_t *nav)
{
prcopt_t *opt=&rtk->opt;
gtime_t time=obs[0].time;
double *rs,*dts,*var,*y,*e,*azel,*freq,*v,*H,*R,*xp,*Pp,*xa,*bias,dt;
int i,j,f,n=nu+nr,ns,ny,nv,sat[MAXSAT],iu[MAXSAT],ir[MAXSAT],niter;
int info,vflg[MAXOBS*NFREQ*2+1],svh[MAXOBS*2];
int stat=rtk->opt.mode<=PMODE_DGPS?SOLQ_DGPS:SOLQ_FLOAT;
int nf=opt->ionoopt==IONOOPT_IFLC?1:opt->nf;
trace(3,"relpos : nx=%d nu=%d nr=%d\n",rtk->nx,nu,nr);
dt=timediff(time,obs[nu].time); // 计算流动站，参考站时间差
rs=mat(6,n); dts=mat(2,n); var=mat(1,n); y=mat(nf*2,n); e=mat(3,n);
azel=zeros(2,n); freq=zeros(nf,n);
for (i=0;i<MAXSAT;i++) {
rtk->ssat[i].sys=satsys(i+1,NULL);
for (j=0;j<NFREQ;j++) rtk->ssat[i].vsat[j]=0;
for (j=1;j<NFREQ;j++) rtk->ssat[i].snr [j]=0;
}
// 根据卫星星历计算当前历元下各卫星的位置、速度、钟差
/* satellite positions/clocks */
satposs(time,obs,n,nav,opt->sateph,rs,dts,var,svh);
// 计算基站的各卫星观测值的非差闭合差（观测值-计算值）及卫星的高度角、方位角、卫星矢量等。
/* UD (undifferenced) residuals for base station */
if (!zdres(1,obs+nu,nr,rs+nu*6,dts+nu*2,var+nu,svh+nu,nav,rtk->rb,opt,1,
y+nu*nf*2,e+nu*3,azel+nu*2,freq+nu*nf)) {
errmsg(rtk,"initial base station position error\n");
free(rs); free(dts); free(var); free(y); free(e); free(azel);
free(freq);
return 0;
}
// 后处理中，需要时，调用 intpres 进行插值
/* time-interpolation of residuals (for post-processing) */
if (opt->intpref) {
dt=intpres(time,obs+nu,nr,nav,rtk,y+nu*nf*2);
}
// 选择基站和流动站之间的同步观测卫星。进行 RTK 算法时只需要基线间的同步观测卫星。
// 返回共同观测的卫星个数，输出卫星号列表 sat、在接收机观测值中的 index 值列表 iu 和在基站观测值中的 index 值列表 ir。
/* select common satellites between rover and base-station */
if ((ns=selsat(obs,azel,nu,nr,opt,sat,iu,ir))<=0) {
errmsg(rtk,"no common satellite\n");
free(rs); free(dts); free(var); free(y); free(e); free(azel);
free(freq);
return 0;
}
// 调用 udstate 更新状态值 rtk->x 及其误差协方差 rtk->P
/* temporal update of states */
udstate(rtk,obs,sat,iu,ir,ns,nav);
// 初始化变量内存以及赋初值
trace(4,"x(0)="); tracemat(4,rtk->x,1,NR(opt),13,4);
xp=mat(rtk->nx,1); Pp=zeros(rtk->nx,rtk->nx); xa=mat(rtk->nx,1);
matcpy(xp,rtk->x,rtk->nx,1);
ny=ns*nf*2+2;
v=mat(ny,1); H=zeros(rtk->nx,ny); R=mat(ny,ny); bias=mat(rtk->nx,1);
// 设置迭代次数 niter，动基线加2次
/* add 2 iterations for baseline-constraint moving-base */
niter=opt->niter+(opt->mode==PMODE_MOVEB&&opt->baseline[0]>0.0?2:0);
// 迭代量测更新 niter 次
for (i=0;i<niter;i++) {
// 计算流动站的各位卫星观测值非差闭合差（观测值-计算值）及卫星的高度角、方位角、卫星矢量等。
/* UD (undifferenced) residuals for rover */
if (!zdres(0,obs,nu,rs,dts,var,svh,nav,xp,opt,0,y,e,azel,freq)) {
errmsg(rtk,"rover initial position error\n");
stat=SOLQ_NONE;
break;
}
// 计算双差残差 v，雅可比矩阵 H、双差误差协方差矩阵 R。
/* DD (double-differenced) residuals and partial derivatives */
if ((nv=ddres(rtk,nav,dt,xp,Pp,sat,y,e,azel,freq,iu,ir,ns,v,H,R,
vflg))<1) {
errmsg(rtk,"no double-differenced residual\n");
stat=SOLQ_NONE;
break;
}
// Kalman滤波量测更新，获得浮点解：
// 在进入 filter 函数后需要对状态向量 x、状态误差协方差矩阵 P、雅可比矩阵 H 进行矩阵简化，去除未使用卫星；
// H 矩阵使用的是常规的定义的雅可比矩阵的转置，这主要是因为 RTKLIB 内存储矩阵的规则是按列存储的。
/* Kalman filter measurement update */
matcpy(Pp,rtk->P,rtk->nx,rtk->nx);
if ((info=filter(xp,Pp,H,v,R,rtk->nx,nv))) {
errmsg(rtk,"filter error (info=%d)\n",info);
stat=SOLQ_NONE;
break;
}
trace(4,"x(%d)=",i+1); tracemat(4,xp,1,NR(opt),13,4);
}
// 量测更新完成，再次调用 zdres 和 ddres 更新双差相位/伪距残差，并调用 valpos 进行验证，
// 若通过则更新 rtk->x 以及 rtk->P，并更新模糊度控制结构体。
if (stat!=SOLQ_NONE&&zdres(0,obs,nu,rs,dts,var,svh,nav,xp,opt,0,y,e,azel,
freq)) {
// 利用浮点结果计算双差残差和量测噪声
/* post-fit residuals for float solution */
nv=ddres(rtk,nav,dt,xp,Pp,sat,y,e,azel,freq,iu,ir,ns,v,NULL,R,vflg);
// 进行浮点解有效性验证
/* validation of float solution */
if (valpos(rtk,v,R,vflg,nv,4.0)) {
// 存储浮点结果
/* update state and covariance matrix */
matcpy(rtk->x,xp,rtk->nx,1);
matcpy(rtk->P,Pp,rtk->nx,rtk->nx);
// 存模糊度相关信息，统计有效卫星数
/* update ambiguity control struct */
rtk->sol.ns=0;
for (i=0;i<ns;i++) for (f=0;f<nf;f++) {
if (!rtk->ssat[sat[i]-1].vsat[f]) continue;
rtk->ssat[sat[i]-1].lock[f]++;
rtk->ssat[sat[i]-1].outc[f]=0;
if (f==0) rtk->sol.ns++; /* valid satellite count by L1 */
}
/* lack of valid satellites */ //检验卫星是否有效
if (rtk->sol.ns<4) stat=SOLQ_NONE;
}
else stat=SOLQ_NONE;
}
// 利用 LAMBDA 算法进行模糊度的固定，得到 RTK 固定解。
/* resolve integer ambiguity by LAMBDA */
if (stat!=SOLQ_NONE&&resamb_LAMBDA(rtk,bias,xa)>1) {
// 模糊度解算成功，根据固定结果计算残差和协方差，并进行校验
if (zdres(0,obs,nu,rs,dts,var,svh,nav,xa,opt,0,y,e,azel,freq)) {
/* post-fit reisiduals for fixed solution */
nv=ddres(rtk,nav,dt,xa,NULL,sat,y,e,azel,freq,iu,ir,ns,v,NULL,R,
vflg);
// 进行固定解有效性验证
/* validation of fixed solution */
if (valpos(rtk,v,R,vflg,nv,4.0)) {
// 固定解验证有效，若为hold模式，需要存模糊度信息
/* hold integer ambiguity */
if (++rtk->nfix>=rtk->opt.minfix&&
rtk->opt.modear==ARMODE_FIXHOLD) {
holdamb(rtk,xa);
}
stat=SOLQ_FIX;
}
}
}
// 保存 solution 状态,位置，速度，方差
/* save solution status */
if (stat==SOLQ_FIX) { // 模糊度固定则记录固定解，固定解在 rtk->xa、rtk->Pa 中。
for (i=0;i<3;i++) {
rtk->sol.rr[i]=rtk->xa[i];
rtk->sol.qr[i]=(float)rtk->Pa[i+i*rtk->na];
}
rtk->sol.qr[3]=(float)rtk->Pa[1];
rtk->sol.qr[4]=(float)rtk->Pa[1+2*rtk->na];
rtk->sol.qr[5]=(float)rtk->Pa[2];
if (rtk->opt.dynamics) { /* velocity and covariance */
for (i=3;i<6;i++) {
rtk->sol.rr[i]=rtk->xa[i];
rtk->sol.qv[i-3]=(float)rtk->Pa[i+i*rtk->na];
}
rtk->sol.qv[3]=(float)rtk->Pa[4+3*rtk->na];
rtk->sol.qv[4]=(float)rtk->Pa[5+4*rtk->na];
rtk->sol.qv[5]=(float)rtk->Pa[5+3*rtk->na];
}
}
else { // 浮点解，浮点解数据在 rtk->x、rtk->P中。
for (i=0;i<3;i++) {
rtk->sol.rr[i]=rtk->x[i];
rtk->sol.qr[i]=(float)rtk->P[i+i*rtk->nx];
}
rtk->sol.qr[3]=(float)rtk->P[1];
rtk->sol.qr[4]=(float)rtk->P[1+2*rtk->nx];
rtk->sol.qr[5]=(float)rtk->P[2];
if (rtk->opt.dynamics) { /* velocity and covariance */
for (i=3;i<6;i++) {
rtk->sol.rr[i]=rtk->x[i];
rtk->sol.qv[i-3]=(float)rtk->P[i+i*rtk->nx];
}
rtk->sol.qv[3]=(float)rtk->P[4+3*rtk->nx];
rtk->sol.qv[4]=(float)rtk->P[5+4*rtk->nx];
rtk->sol.qv[5]=(float)rtk->P[5+3*rtk->nx];
}
rtk->nfix=0;
}
// 存当前历元载波信息，供下次使用
for (i=0;i<n;i++) for (j=0;j<nf;j++) {
if (obs[i].L[j]==0.0) continue;
rtk->ssat[obs[i].sat-1].pt[obs[i].rcv-1][j]=obs[i].time; //先前历元载波相位时间
rtk->ssat[obs[i].sat-1].ph[obs[i].rcv-1][j]=obs[i].L[j]; //先前历元载波相位观测值
}
// 存 SNR 信噪比信息
for (i=0;i<ns;i++) for (j=0;j<nf;j++) {
/* output snr of rover receiver */
rtk->ssat[sat[i]-1].snr[j]=obs[iu[i]].SNR[j];
}
// 存卫星的 fix 信息及周跳信息
for (i=0;i<MAXSAT;i++) for (j=0;j<nf;j++) {
if (rtk->ssat[i].fix[j]==2&&stat!=SOLQ_FIX) rtk->ssat[i].fix[j]=1;
if (rtk->ssat[i].slip[j]&1) rtk->ssat[i].slipc[j]++;
}
free(rs); free(dts); free(var); free(y); free(e); free(azel); free(freq);
free(xp); free(Pp); free(xa); free(v); free(H); free(R); free(bias);
if (stat!=SOLQ_NONE) rtk->sol.stat=stat;
return stat!=SOLQ_NONE;
}9.2.3 intpres()：基站残差的插值（后处理）
intpres() 是一个不重要的函数，它仅用于后处理。
1. 参数列表
c
/* args */
gtime_t time I 历元时间
const obsd_t *obs I 卫星原始观测数据
int n I 观测数据数目
const nav_t *nav I 星历数据
rtk_t *rtk IO rtk 控制/状态结构体
double *y IO 基站非差残差
/* return */
int time - fabs(ttb)<fabs(tt)?ttb:tt2. 执行流程
- 时间差计算：
tt = timediff(time, obs[0].time)：计算流动站时间与当前基站观测时间的差值。ttb = timediff(time, obsb[0].time)：计算流动站时间与之前基站观测时间的差值。
- 卫星位置计算：
- 调用
satposs()函数，基于之前基站观测时间obsb计算卫星位置（rs）、钟差（dts）、方差（var）和健康状态（svh）。
- 调用
- 残差计算：
- 调用
zdres()函数，计算之前基站观测的伪距减去几何距离的残差（yb），并获取卫星仰角/方位角（azel）和频率（freq）。 - 若
zdres()失败，返回tt。
- 调用
- 残差插值：
- 遍历当前基站观测的每个卫星（
i从 0 到 n-1）。 - 找到对应之前基站观测中的卫星（j），若无匹配则跳过。
- 对于每个频率（
k从 0 到 nf*2-1）：p指向当前残差y[i*nf*2+k]，q指向之前残差yb[j*nf*2+k]。- 若当前或之前残差为 0，或存在周跳（LLI & LLI_SLIP），则将当前残差置为 0。
- 否则，使用线性插值公式
*p = (ttb*(*p) - tt*(*q))/(ttb-tt)计算插值后的残差。
- 遍历当前基站观测的每个卫星（
- 返回时间差：
- 返回绝对值较小的时间差
（fabs(ttb)<fabs(tt)?ttb:tt），表示插值参考时间。
- 返回绝对值较小的时间差
3. 注意事项
- 由于观测方程的本质是最小二乘，而最小二乘中的残差项包含了观测信息。所以我们可以说残差的插值，就是对观测数据的插值，只是计算顺序上就是先处理观测信息，再进行 OMC（Observation Minus Computed）运算，但二者本质上是等价的。
4. 源码注释
点击查看代码
c
/* time-interpolation of residuals (for post-processing solutions) -----------
time = rover time stamp
obs = pointer to first base observation for this epoch
y = pointer to base obs errors */
static double intpres(gtime_t time, const obsd_t *obs, int n, const nav_t *nav,
rtk_t *rtk, double *y)
{
static obsd_t obsb[MAXOBS]; /* 存储之前基站观测数据的静态数组 */
static double yb[MAXOBS*NFREQ*2], rs[MAXOBS*6], dts[MAXOBS*2], var[MAXOBS]; /* 之前残差、卫星位置、钟差、方差 */
static double e[MAXOBS*3], azel[MAXOBS*2], freq[MAXOBS*NFREQ]; /* 误差、仰角/方位角、频率 */
static int nb = 0, svh[MAXOBS*2]; /* 之前基站卫星数、卫星健康状态 */
prcopt_t *opt = &rtk->opt; /* 解算选项指针 */
double tt, ttb, *p, *q; /* 时间差、指针 */
int i, j, k, nf = NF(opt); /* 循环索引、频率数 */
tt = timediff(time, obs[0].time); /* 计算流动站与当前基站观测时间差 */
trace(3, "intpres : n=%d tt=%.1f, epoch=%d\n", n, tt, rtk->epoch); /* 调试日志：输出卫星数、时间差、历元 */
/* 若为第一个历元或时间差很小，直接使用当前基站观测更新之前数据 */
if (nb == 0 || rtk->epoch == 0 || fabs(tt) < DTTOL) {
nb = n; /* 更新之前卫星数 */
for (i = 0; i < n; i++) obsb[i] = obs[i]; /* 当前基站观测复制到之前数组 */
return tt; /* 返回时间差 */
}
/* 若与之前基站时间差过大或与当前时间差相同，不进行插值 */
ttb = timediff(time, obsb[0].time); /* 计算流动站与之前基站观测时间差 */
if (fabs(ttb) > opt->maxtdiff * 2.0 || ttb == tt) return tt;
/* 计算之前基站观测的卫星位置 */
satposs(time, obsb, nb, nav, opt->sateph, rs, dts, var, svh);
/* 计算之前基站观测的伪距减几何距离的残差 */
if (!zdres(1, obsb, nb, rs, dts, var, svh, nav, rtk->rb, opt, yb, e, azel, freq)) {
return tt; /* 若计算失败，返回时间差 */
}
/* 插值当前和之前基站观测的残差 */
for (i = 0; i < n; i++) {
/* 匹配当前卫星与之前卫星 */
for (j = 0; j < nb; j++) if (obsb[j].sat == obs[i].sat) break;
if (j >= nb) continue; /* 若无匹配，跳过 */
/* p 指向当前残差，q 指向之前残差 */
for (k = 0, p = y + i * nf * 2, q = yb + j * nf * 2; k < nf * 2; k++, p++, q++) {
if (*p == 0.0 || *q == 0.0 || /* 若残差为 0 */
(obs[i].LLI[k % nf] & LLI_SLIP) || /* 当前观测有周跳 */
(obsb[j].LLI[k % nf] & LLI_SLIP)) { /* 之前观测有周跳 */
*p = 0.0; /* 置当前残差为 0 */
} else {
/* 线性插值：(ttb*当前 - tt*之前)/(ttb-tt) */
*p = (ttb * (*p) - tt * (*q)) / (ttb - tt);
}
}
}
/* 返回较小的时间差作为插值参考 */
return fabs(ttb) < fabs(tt) ? ttb : tt;
}9.2.4 selsat()：选择共视卫星
selsat() 函数用于共视星的选取，它如何实现反倒不重要。
1. 参数列表
c
/* args */
const obsd_t *obs I 卫星原始观测数据
double *azel I 卫星 [方位角，仰角]
int nu I 接收机观测数据数目
int nr I 基站观测数据数目
const prcopt_t *opt I 配置选项
int *sat O 共视卫星数组
int *iu O 流动站共视卫星索引
int *ir O 基站的公式卫星索引
/* return */
int ns - 共视卫星的个数2. 注意事项
- 返回基站和流动站之间的共视星个数
ns，输出卫星号列表sat，及流动站的索引列表iu和在基站的索引列表ir obs[i]中：流动站为obs[0~nu-1]，基站为obs[nu~nu+nr-1]
3. 源码注释
点击查看代码
c
static int selsat(const obsd_t *obs, double *azel, int nu, int nr,
const prcopt_t *opt, int *sat, int *iu, int *ir)
{
int i,j,k=0;
trace(3,"selsat : nu=%d nr=%d\n",nu,nr);
for (i=0,j=nu;i<nu&&j<nu+nr;i++,j++) {
if (obs[i].sat<obs[j].sat) j--;
else if (obs[i].sat>obs[j].sat) i--;
else if (azel[1+j*2]>=opt->elmin) { /* elevation at base station */
sat[k]=obs[i].sat; iu[k]=i; ir[k++]=j;
trace(4,"(%2d) sat=%3d iu=%2d ir=%2d\n",k-1,obs[i].sat,i,j);
}
}
return k;
}9.2.5 valpos()：解的有效性验证
计算流动站非差、双差残差，残差平方和 v[i]*v[i] 是否小于 fact*R[i+i*nv]。
1. 参数列表
c
/* args */
rtk_t *rtk IO 解结构体
const double *v I 残差
const double *R I 观测误差协方差
const int *vflg I 双差中使用的卫星列表
int nv I 残差的个数
double thres I 残差阈值
/* return */
int stat - 状态（1：有效 0：无效）2. 注意事项
- 以 demo5 的代码为基础，无论是单点定位还是相对定位，
valpos()函数更多的只是用来评估和输出解的状态信息，并没有实际影响代码的执行。它可能会检测出残差过大的卫星，但不会对卫星做任何后续处理。
3. 源码注释
点击查看代码
c
static int valpos(rtk_t *rtk, const double *v, const double *R, const int *vflg,
int nv, double thres)
{
double fact=thres*thres;
int i,stat=1,sat1,sat2,type,freq;
char *stype;
trace(3,"valpos : nv=%d thres=%.1f\n",nv,thres);
/* post-fit residual test */
for (i=0;i<nv;i++) {
if (v[i]*v[i]<=fact*R[i+i*nv]) continue;
sat1=(vflg[i]>>16)&0xFF; // 参考卫星号
sat2=(vflg[i]>> 8)&0xFF; // 非参考卫星号
type=(vflg[i]>> 4)&0xF; // 种类：0：载波；1：伪距；3：动基线
freq=vflg[i]&0xF; // 第几个频率（从0开始）
stype=type==0?"L":(type==1?"L":"C");
errmsg(rtk,"large residual (sat=%2d-%2d %s%d v=%6.3f sig=%.3f)\n",
sat1,sat2,stype,freq+1,v[i],SQRT(R[i+i*nv]));
}
return stat;
}9.2.6 outsolstat()：输出解状态信息
解状态信息除了参考其内部调用的 rtkoutstat() 函数，也可以关注相应的格式定义，具体参考 RTKLIB-Manual-CN B.3节：解状态文件。
1. 参数列表
c
/* args */
rtk_t *rtk IO rtk 控制/结果结构体
const nav_t *nav I 卫星星历数据
/* return */
None2. 内容明细
a. 位置状态参数
$POS,week,tow,stat,posx,posy,posz,posxf,posyf,poszfweek/tow: GPS 周和周内秒（s）stat: 定位结果状态posx/posy/posz: 浮点解坐标 x/y/z ecef (m)posxf/posyf/poszf: 固定解坐标 x/y/z ecef (m)
b. 速度、加速度状态参数
$VELACC,week,tow,stat,vele,veln,velu,acce,accn,accu,velef,velnf,veluf,accef,accnf,accufweek/tow: GPS 周和周内秒（s）stat: 定位结果状态vele/veln/velu: 浮点解速度 e/n/u (m/s)acce/accn/accu: 浮点解加速度 e/n/u (m/s^2)velef/velnf/veluf: 固定解速度 e/n/u (m/s)accef/accnf/accuf: 固定解加速度 e/n/u (m/s^2)
c. 接收机钟差状态参数
$CLK,week,tow,stat,rcv,clk1,clk2,clk3,clk4week/tow: GPS 周和周内秒（s）stat: 定位结果状态rcv: 接收机 (1:rover,2:base station)clk1: 接收机 GPS 钟差 (ns)clk2: 接收机 GLONASS 钟差 (ns)clk3: 保留clk4: 保留
d. 估计电离层状态参数
$ION,week,tow,stat,sat,az,el,ion,ion‐fixedweek/tow: GPS 周和周内秒（s）stat: 定位结果状态sat: 卫星 idaz/el: 方位角/仰角 (deg)ion: 浮点解垂直电离层延迟 L1 (m)ion‐fixed: 固定解垂直电离层延迟 L1 (m)
e. 估计对流层状态参数
$TROP,week,tow,stat,rcv,ztd,ztdfweek/tow: GPS 周和周内秒（s）stat: 定位结果状态rcv: 接收机标识 (1:rover,2:base station)ztd: 浮点解总对流延迟（m）ztdf: 固定解总对流延迟（m）
f. 估计GLONASS receiver H/W bias difference参数
$HWBIAS,week,tow,stat,frq,bias,biasfweek/tow: GPS 周和周内秒（s）stat: 定位结果状态frq: 频率（1:L1,2:L2,3:L5,...）bias: 接收机 H/W 偏差系数（m/MHz）浮点解biasf: 接收机 H/W 偏差系数（m/MHz）固定解
g. 伪距和载波相位观测量的残差
$SAT,week,tow,sat,frq,az,el,resp,resc,vsat,snr,fix,slip,lock,outc,slipc,rejcweek/tow: GPS 周和周内秒（s）sat/frq: 卫星 id/频率（1:L1,2:L2,3:L5,...）az/el: 方位角/仰角（deg）resp: 伪距残差（m）resc: 载波相位残差（m）vsat: 有效数据标志（0:无效,1:有效）snr: 信号强度（dbHz）fix: 模糊度固定状态标识 (0:no data,1:float,2:fixed,3:hold)slip: 周跳标识 (bit1:slip,bit2:parity unknown)lock: 载波观测锁定计数outc: 数据失效标识计数slipc: 周跳计数rejc: 数据拒绝计数（异常值）
3. 源码注释
点击查看代码
c
static void outsolstat(rtk_t *rtk)
{
ssat_t *ssat;
double tow;
char buff[MAXSOLMSG+1],id[32];
int i,j,n,week,nfreq,nf=NF(&rtk->opt);
if (statlevel<=0||!fp_stat||!rtk->sol.stat) return;
trace(3,"outsolstat:\n");
// 根据时间分结果文件
/* swap solution status file */
swapsolstat();
/* write solution status */
n=rtkoutstat(rtk,buff); buff[n]='\0';
fputs(buff,fp_stat);
// 如果解的状态为 SOLQ_NONE，或结果输出等级小于等于 1，直接直接返回
if (rtk->sol.stat==SOLQ_NONE||statlevel<=1) return;
tow=time2gpst(rtk->sol.time,&week);
nfreq=rtk->opt.mode>=PMODE_DGPS?nf:1;
/* write residuals and status */
for (i=0;i<MAXSAT;i++) {
ssat=rtk->ssat+i;
if (!ssat->vs) continue;
satno2id(i+1,id);
for (j=0;j<nfreq;j++) {
fprintf(fp_stat,"$SAT,%d,%.3f,%s,%d,%.1f,%.1f,%.4f,%.4f,%d,%.1f,%d,%d,%d,%d,%d,%d\n",
week,tow,id,j+1,ssat->azel[0]*R2D,ssat->azel[1]*R2D,
ssat->resp[j],ssat->resc[j],ssat->vsat[j],
ssat->snr[j]*SNR_UNIT,ssat->fix[j],ssat->slip[j]&3,
ssat->lock[j],ssat->outc[j],ssat->slipc[j],ssat->rejc[j]);
}
}
}