对于pppos中Kalman滤波的思考
假设:
H 26*216 几何矩阵
P 216*216 参数预测值的协方差阵
R 26*26 观测噪声的协方差阵
x 216*1 参数预测值
xp 216*1 参数平差值
Pp 216*216 参数平差值得协方差阵
v 26*1 新息向量
增益矩阵
静态的情况下,RTKlib还是进行的是Kalman滤波,并不是最小二乘。
因为对于Kalman滤波,其通用解的形式是:
是参数预测值的协方差阵,并不是为0。
对于pos
在上一个历元的参数平差值的协方差阵+过程噪声(过程噪声为0)
/*位置参数也是上一个历元的平差值*/
if (rtk->opt.mode==PMODE_PPP_STATIC) {
for (i=0;i<3;i++) {
rtk->P[i*(1+rtk->nx)]+=SQR(rtk->opt.prn[5])*fabs(rtk->tt);
}
return;}
对于clk
/*计算钟差的协方差阵,更新到rtk.P中*/
/*钟差是上一个历元的平差值*/
initx(rtk,CLIGHT*dtr,VAR_CLK,IC(i,&rtk->opt));
对于对流层:在上一个历元的参数平差值的协方差阵+过程噪声
rtk->P[i+i*rtk->nx]+=SQR(rtk->opt.prn[2])*fabs(rtk->tt);
对于模糊度参数:
/*offset/k 平均值*/每颗卫星的bias与上一个历元平差后的模糊度参数rtk->x之间的偏差offset之和
if (k>=2&&fabs(offset/k)>0.0005*CLIGHT) {
for (i=0;i<MAXSAT;i++) {
j=IB(i+1,f,&rtk->opt);
if (rtk->x[j]!=0.0) rtk->x[j]+=offset/k;改变了,上一历元的平差值+此历元计算出来的偏移之和平均值
}
trace(2,"phase-code jump corrected: %s n=%2d dt=%12.9fs\n",
time_str(rtk->sol.time,0),k,offset/k/CLIGHT);
}
for (i=0;i<n&&i<MAXOBS;i++) {
sat=obs[i].sat;
j=IB(sat,f,&rtk->opt);
/*prn[0] bias process-noise std*/
rtk->P[j+j*rtk->nx]+=SQR(rtk->opt.prn[0])*fabs(rtk->tt); 改变了,上一历元的平差值的协方差阵+过程噪声(为0)
注意:这里的bias[i]==0.0 是不是那些对于无电离层组合不敏感的周跳
if (bias[i]==0.0||(rtk->x[j]!=0.0&&!slip[i])) continue;
/* reinitialize phase-bias if detecting cycle slip */
initx(rtk,bias[i],VAR_BIAS,IB(sat,f,&rtk->opt));
/* reset fix flags */
for (k=0;k<MAXSAT;k++) rtk->ambc[sat-1].flags[k]=0;
trace(5,"udbias_ppp: sat=%2d bias=%.3f\n",sat,bias[i]);
}/*end sat */