Winchell

Appearance

9.3 udstate(): 时间更新 ​

udrcvbias() 用于 GLONASS 接收机进行硬件偏移的时间更新操作，这部分为 demo5 新增内容，并非 RTKLIB 的重点，因此暂不进行解析。

9.3.1 udstate(): 函数主体 ​

图9.2-2 udstate() 函数调用

Kalman 滤波的时间更新，更新状态值 rtk->x 及其误差协方差 rtk->P。

1. 参数列表

c
/* args */
rtk_t    *rtk      IO  rtk控制结构体
obsd_t   *obs      I   观测数据
int       sat      I   接收机和基站共同观测的卫星号列表
int      *iu       I   移动站共视星所在obs数组中的索引
int      *ir       I   基站共视星所在obs数组中的索引
int       ns       I   接收机和基站共同观测的卫星个数
nav_t    *nav      I   导航数据
/* return */
none

2. 执行流程

  • PVA 更新：调用 udpos() 根据不同模式更新 rtk 中的位置、速度、加速度值和协方差。
  • 电离层参数更新：电离层模式 >=IONOOPT_EST，调用 baseline() 计算基线长度 BL，调用 udion() 根据基线长度 bl 更新状态 rtk->x 中的电离层参数（MAXSAT 个）及其协方差。

    在基线比较长的情况下（>10km）,电离层的误差无法被双差完全消除，因此需要考虑电离层的影响。此时可以在配置中选择电离层修正方式为 IONOOPT_EST。这种配置下，会将垂直方向的单差电离层延迟添加到卡尔曼滤波状态量中。

  • 对流层参数更新：若对流层模式 >=TROPOPT_EST，调用 udtrop() 根据基线长度 BL 更新状态 rtk->x 中的对流层参数（2或6个）及其协方差。

    调用 udtrop() 会传入电离层模式 >=IONOOPT_EST，调用 baseline() 计算基线长度 BL。

  • 接收机硬件偏移更新：若为 GLONASS AR 模式，调用 udrcvbias() 更新接收机硬件偏移（demo5代码新增，仅针对 u-blox 接收机）。
  • 单差模糊度更新：若模式 >PMODE_DGPS，则使用载波相位定位，调用 udbias() 更新载波相位偏移状态值以及其误差协方差（udbias 内部还会进行周跳探测，这部分参考9.7节）。

3. 源码注释

点击查看代码
c
static void udstate(rtk_t *rtk, const obsd_t *obs, const int *sat,
                    const int *iu, const int *ir, int ns, const nav_t *nav)
{
    double tt=rtk->tt,bl,dr[3];
    
    trace(3,"udstate : ns=%d\n",ns);

    // 调用 udpos 根据不同模式更新 rtk 中的位置、速度、加速度值和协方差
    /* temporal update of position/velocity/acceleration */
    udpos(rtk,tt);
    
    // 电离层模式 >=IONOOPT_EST，调用 udion 更新状态 rtk->x 中的电离层参数（MAXSAT 个）及其协方差
    /* temporal update of ionospheric parameters */
    if (rtk->opt.ionoopt>=IONOOPT_EST) {    
        bl=baseline(rtk->x,rtk->rb,dr);
        udion(rtk,tt,bl,sat,ns);
    }

    // 若对流层模式 >=TROPOPT_EST，调用 udtrop 更新状态 rtk->x 中的对流层参数（2或6个）及其协方差
    /* temporal update of tropospheric parameters */
    if (rtk->opt.tropopt>=TROPOPT_EST) {
        udtrop(rtk,tt,bl); // 也需要用到上面电离层处理时的基线长度bl
    }

    // 若为 GLONASS AR 模式，调用 udrcvbias 更新接收机硬件偏移。
    /* temporal update of eceiver h/w bias */
    if (rtk->opt.glomodear==2&&(rtk->opt.navsys&SYS_GLO)) {
        udrcvbias(rtk,tt);
    }

    // 若模式 >PMODE_DGPS，调用 udbias 更新载波相位偏移状态值以及其误差协方差。
    /* temporal update of phase-bias */
    if (rtk->opt.mode>PMODE_DGPS) {
        udbias(rtk,tt,obs,sat,iu,ir,ns,nav);
    }
}

9.3.2 baseline()：求基线长度，流动站坐标减基站坐标的模 ​

1. 参数列表

c
/* args */
const double  *ru    I   流动站坐标
const double  *rb    I   基站坐标
double        *dr    IO  基线矢量的单位向量
/* return */
double         bl    O   基线长度

2. 执行流程

  • 计算基线矢量 dr，返回基线长度 bl。
c
static double baseline(const double *ru, const double *rb, double *dr)
{
    int i;
    for (i=0;i<3;i++) dr[i]=ru[i]-rb[i];
    return norm(dr,3);
}

9.3.3 udpos()：位置参数时间更新 ​

更新 rtk 中的位置、速度、加速度值和协方差。

1. 流动站运动模式：

  • 固定模式（已知流动站坐标）：以流动站坐标为每一历元的时间更新值，并给协方差阵 P 设置一个较小的方差值 (1E-8)。
  • 静态定位：以上一历元的解作为这一历元的时间更新值，不进行更新操作。
  • 动态定位（非动态模式）：以这一历元的单点定位解作为位置状态量的时间更新值，方差为默认 (30*30)，状态中不包含速度和加速度。
  • 动态定位（动态模式）：假设历元间保持匀速运动，根据速度、加速度完成时间更新，协方差过大会用伪距、多普勒计算的位置、速度重置，加速度设为极小值 1E-6。

2. 参数列表：

c
/* args */
rtk_t    *rtk      IO    rtk 控制/解算结果结构体
double    tt       I     本次更新与上次更新的时间差
/* return */
none

3. 执行流程：

  • PMODE_FIXED 模式的更新：若为 PMODE_FIXED 模式，直接从选项中取得位置值赋值给 rtk->x，并将其协方差阵 P 的方差设置为一个较小的方差值，然后直接返回；
  • 第一个历元的初始化：若为第一个历元，用 rtk->sol 中的位置值初始化 rtk->x，位置状态量初始化为单点定位所得到的位置值，方差设置为 VAR_POS（认为初值误差较大，需要迭代收敛）；
  • PMODE_STATIC 模式的更新：若为 PMODE_STATIC 静态模式，上一时刻状态就是当前时刻状态；
  • 非动态模式的更新（dynamics = off，状态中只有位置）：
    • 若不是 dynamics 模式，那么每次状态更新都将位置状态量赋值为单点定位的位置，方差设置为 VAR_POS。
  • 动态模式的更新（dynamics = on，状态中有位置、速度和加速度）：

    • 位置方差的检验：如果位置状态量的平均方差大于阈值 VAR_POS 则用 rtk->sol 中的位置值、速度值重置 rtk->x[0~2]、rtk->x[3~5]，1E-6赋值 rtk->x[3~5]，方差分别设为 VAR_POS、VAR_VEL、VAR_ACC；

    • 矩阵去除非零元素：生成有效状态（非零元素）下标 ix[]；

    • F 阵的构建：用时间差 tt 构建状态转移矩阵 F；

      其中：

      这里的 对应了源码中的 F，如果考虑了接收机动态（pos1-dynamics=on），那么还需要包含加速度（常加速度模型）。

      上面的公式描述了 EKF 系统状态递推的过程，其中的 ，表示采样时间间隔。

    • rtk->x 赋值参数阵 x[]、rtk->P 赋值方差阵 P[]；

    • Kalman 时间更新：状态更新 rtk->x，rtk->P，, ；

    • Q 阵的构建及 P 阵的更新：给加速度加过程噪声(随机游走噪声) rtk->opt.prn[3]，rtk->opt.prn[4]，先存到 Q[]，再赋值给 rtk->P。

4. 注意事项

  • initx() 的使用：赋值 xi、var、将 rtk->x[i]，rtk->P[i,i]，rtk->P[i] 上非对角线元素赋值为 0。虽然函数中有初始化 init 字样，但是其实 initx() 会在很多赋值的场合下使用。

5. 源码注释

点击查看代码
c
static void udpos(rtk_t *rtk, double tt)
{
    double *F,*P,*FP,*x,*xp,pos[3],Q[9]={0},Qv[9],var=0.0;
    int i,j,*ix,nx;
    
    trace(3,"udpos   : tt=%.3f\n",tt);
    // 若为 PMODE_FIXED 模式，直接从选项中取得位置值给 rtk->x，然后返回
    /* fixed mode */
    if (rtk->opt.mode==PMODE_FIXED) {
        for (i=0;i<3;i++) initx(rtk,rtk->opt.ru[i],1E-8,i);
        return;
    }
    
    // 若为第一个历元，用 rtk->sol 中的位置值初始化 rtk->x
    /* 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);
        if (rtk->opt.dynamics) {
            for (i=3;i<6;i++) initx(rtk,rtk->sol.rr[i],VAR_VEL,i);
            for (i=6;i<9;i++) initx(rtk,1E-6,VAR_ACC,i);
        }
    }
    /* static mode */
    if (rtk->opt.mode==PMODE_STATIC) return; // 若为 PMODE_STATIC 模式，return
    
    // 若非 dynamics 模式，用 rtk->sol 中的位置值初始化 rtk->x，return
    /* kinmatic mode without dynamics */
    if (!rtk->opt.dynamics) {
        for (i=0;i<3;i++) initx(rtk,rtk->sol.rr[i],VAR_POS,i);
        return;
    }

    // 检查位置协方差，若大于阈值 VAR_POS 则用 rtk->sol 中的位置值重置 rtk->x
    /* check variance of estimated postion */
    for (i=0;i<3;i++) var+=rtk->P[i+i*rtk->nx];
    var/=3.0;
    if (var>VAR_POS) {
        /* reset position with large variance */
        for (i=0;i<3;i++) initx(rtk,rtk->sol.rr[i],VAR_POS,i);  //位置
        for (i=3;i<6;i++) initx(rtk,rtk->sol.rr[i],VAR_VEL,i);  //速度
        for (i=6;i<9;i++) initx(rtk,1E-6,VAR_ACC,i);            //加速度
        trace(2,"reset rtk position due to large variance: var=%.3f\n",var);
        return;
    }

    // 生成有效状态(非0)下标ix[]
    /* generate valid state index */
    ix=imat(rtk->nx,1);
    for (i=nx=0;i<rtk->nx;i++) {
        if (rtk->x[i]!=0.0&&rtk->P[i+i*rtk->nx]>0.0) ix[nx++]=i;
    }
    if (nx<9) {
        free(ix);
        return;
    }
    
    // 用时间差 tt 构建状态转移矩阵 F
    /* state transition of position/velocity/acceleration */
    F=eye(nx); P=mat(nx,nx); FP=mat(nx,nx); x=mat(nx,1); xp=mat(nx,1);
    for (i=0;i<6;i++) {
        F[i+(i+3)*nx]=tt;
    }
    for (i=0;i<3;i++) {
        F[i+(i+6)*nx]=SQR(tt)/2.0;
    }

    // rtk->x赋值参数阵 x[]、rtk->P 赋值方差阵 P[]
    for (i=0;i<nx;i++) {
        x[i]=rtk->x[ix[i]];
        for (j=0;j<nx;j++) {
            P[i+j*nx]=rtk->P[ix[i]+ix[j]*rtk->nx];
        }
    }

    // 状态更新
    /* x=F*x, P=F*P*F+Q */
    matmul("NN",nx,1,nx,1.0,F,x,0.0,xp);
    matmul("NN",nx,nx,nx,1.0,F,P,0.0,FP);
    matmul("NT",nx,nx,nx,1.0,FP,F,0.0,P);
    for (i=0;i<nx;i++) {
        rtk->x[ix[i]]=xp[i];
        for (j=0;j<nx;j++) {
            [ix[i]+ix[j]*rtk->nx]=P[i+j*nx];
        }
    }
    
    // 给加速度加过程噪声,先存到 Q，再赋值给 rtk->P
    /* process noise added to only acceleration */
    Q[0]=Q[4]=SQR(rtk->opt.prn[3])*fabs(tt);
    Q[8]=SQR(rtk->opt.prn[4])*fabs(tt);
    ecef2pos(rtk->x,pos);
    covecef(pos,Q,Qv);
    for (i=0;i<3;i++) for (j=0;j<3;j++) {
        rtk->P[i+6+(j+6)*rtk->nx]+=Qv[i+j*3];
    }
    free(ix); free(F); free(P); free(FP); free(x); free(xp);
}

9.3.4 udion()：电离层参数时间更新 ​

1. 参数列表：

c
/* args */
IO  rtk_t   *rtk:     rtk solution structure
I   double   tt:      当前历元与前一历元的时间差
I   double   bl:      基线长度
I   int     *sat:     移动站、基站共视星列表
I   int      ns:      共视星个数
/* return */
none

2. 执行流程：

  • 中断超限异常处理：遍历每颗卫星，如果两个频率，载波中断计数都大于GAP_RESION（120），电离层状态量（单差电离层延迟状态量）设为 0；
  • 电离层参数更新：
    • 如果电离层状态量为 0，状态设为 1E-6，协方差设为 SQR(rtk->opt.std[1]*bl/1E4)，其中 bl 为基线长度；
    • 电离层状态量不为 0，加过程噪声 SQR(rtk->opt.prn[1]*bl/1E4*fact)*fabs(tt)。

3. 注意事项

  • IONOOPT_EST 模式才会执行此函数，电离层参数数为卫星数；
  • 长基线情形下，电离层误差无法被双差完全消除，可以配置 IONOOPT_EST 将垂直方向的单差电离层延迟添加到卡尔曼滤波状态量中；
  • 在整个电离层状态量的时间更新中，实际上仅对电离层状态量为 0 的那些卫星进行了初始化，其余卫星的电离层状态量并没有变化，仅仅是在协方差阵中加入了过程噪声。

4. 源码注释

点击查看代码
c
static void udion(rtk_t *rtk, double tt, double bl, const int *sat, int ns)
{
    double el,fact;
    int i,j;
    
    trace(3,"udion   : tt=%.3f bl=%.0f ns=%d\n",tt,bl,ns);
    
    //遍历每颗卫星
    for (i=1;i<=MAXSAT;i++) {
        j=II(i,&rtk->opt);  //获取电离层参数下标j
        //如果两个频率，载波中断计数都大于GAP_RESION(120)，电离层参数设为0
        if (rtk->x[j]!=0.0&&
            rtk->ssat[i-1].outc[0]>GAP_RESION&&rtk->ssat[i-1].outc[1]>GAP_RESION)
            rtk->x[j]=0.0;
    }
    //遍历共视卫星
    for (i=0;i<ns;i++) {
        j=II(sat[i],&rtk->opt); //获取电离层参数下标
        
        //如果电离层状态量为0，状态设为1E-6，协方差设为SQR(rtk->opt.std[1]*bl/1E4)，bl为基线长度
        if (rtk->x[j]==0.0) {
            initx(rtk,1E-6,SQR(rtk->opt.std[1]*bl/1E4),j);
        }
        //电离层状态量不为0，加过程噪声SQR(rtk->opt.prn[1]*bl/1E4*fact)*fabs(tt)
        else {
            /* elevation dependent factor of process noise */
            el=rtk->ssat[sat[i]-1].azel[1];
            fact=cos(el);
            rtk->P[j+j*rtk->nx]+=SQR(rtk->opt.prn[1]*bl/1E4*fact)*fabs(tt);
        }
    }
}

9.3.5 udtrop()：对流层参数时间更新 ​

1. 参数列表：

c
/* args */
rtk_t   *rtk     IO     rtk 控制/解算结果结构体
double   tt      I      当前历元与前一历元的时间差
double   bl      I      基线长度
/* return */
none

2. 执行流程：

  • for 循环，两次分别处理基站和移动站：
    • 如果对流层状态量为0：用INIT_ZWD、SQR(rtk->opt.std[2])初始化天顶方向对流层延迟，TROPOPT_ESTG模式，用1E-6、VAR_GRA初始化东向、北向的对流层梯度系数；
    • 如果对流层状态量不为0：增加过程噪声SQR(rtk->opt.prn[2])*fabs(tt)，TROPOPT_ESTG模式，给东向、北向的对流层梯度系数加过程噪声SQR(rtk->opt.prn[2]*0.3)*fabs(tt)。

3. 注意事项

  • TROPOPT_ESTG 和 TROPOPT_EST 模式才执行此函数。对流层参数数与卫星数无关，TROPOPT_EST 模式为 2 个，TROPOPT_ESTG 模式为 6 个；
  • 在整个对流层状态量的时间更新中，实际上仅对对流层状态量为 0 的那些卫星进行了初始化，其余卫星的对流层状态量并没有变化，仅仅是在协方差阵中加入了过程噪声。
  • 在基线比较长（>10km）、或者是基线和移动站之间高度差较大的情况下，由于对流层误差无法被双差完全消除，因此需要考虑对流层延迟的影响。此时可以在配置中选择对流层修正方式为 Estiamte ZTD 或者 Estimate ZTD+Grad 模式：
    • Estiamte ZTD：将基站、移动站天顶方向的对流层延迟(,)加入卡尔曼滤波状态量，即新增2个状态量；
    • Estimate ZTD+Grad：除了基线、移动天顶方向对流层延迟外，还会将东向、北向的对流层梯度系数,,,加入卡尔曼滤波状态量，即新增6个状态量。

4. 源码注释

点击查看代码
c
static void udtrop(rtk_t *rtk, double tt, double bl)
{
    int i,j,k;
    
    trace(3,"udtrop  : tt=%.3f\n",tt);
    
    //处理基站、移动站对流层延迟
    for (i=0;i<2;i++) { 
        j=IT(i,&rtk->opt);  //获取对流层参数下标j
        //如果对流层状态量为0，
        if (rtk->x[j]==0.0) {   
        //用INIT_ZWD、SQR(rtk->opt.std[2])初始化天顶方向对流层延迟
            initx(rtk,INIT_ZWD,SQR(rtk->opt.std[2]),j); 
            //TROPOPT_ESTG模式，用1E-6、VAR_GRA初始化东向、北向的对流层梯度系数
            if (rtk->opt.tropopt>=TROPOPT_ESTG) {   
                for (k=0;k<2;k++) initx(rtk,1E-6,VAR_GRA,++j);
            }
        }
        //对流层状态量不为0
        else {
            //增加过程噪声SQR(rtk->opt.prn[2])*fabs(tt)
            rtk->P[j+j*rtk->nx]+=SQR(rtk->opt.prn[2])*fabs(tt);
            //TROPOPT_ESTG模式，给东向、北向的对流层梯度系数加过程噪声`SQR(rtk->opt.prn[2]*0.3)*fabs(tt)`
            if (rtk->opt.tropopt>=TROPOPT_ESTG) {   
                for (k=0;k<2;k++) {
                    rtk->P[++j*(1+rtk->nx)]+=SQR(rtk->opt.prn[2]*0.3)*fabs(tt);
                }
            }
        }
    }
}

9.3.6 udbias()：单差相位偏移时间更新（单差整周模糊度） ​

1. 基本原理

udbias 利用单差伪距和单差载波相位计算一个单差相位偏移平均值，从而实现单差模糊度的状态更新，整个过程可以概括为：

偏移估计：

平均偏移：

状态更新：

2. 参数列表

c
/* args */
rtk_t    *rtk      IO  rtk 控制结构体
double    tt       I   本次更新与上次更新的时间差
obsd_t   *obs      I   观测数据
int       sat      I   接收机和基站共同观测的卫星号列表
int      *iu       I   接收机和基站共同观测的卫星在接收机观测值中的 index 值列表
int      *ir       I   接收机和基站共同观测的卫星在基站观测值中的 index 值列表
int       ns       I   接收机和基站共同观测的卫星个数
nav_t    *nav      I   导航数据
/* return */
none

3. 执行流程

  • 周跳探测：首先是对所有共视星进行周跳检测，遍历共视卫星（有关周跳探测更详细的内容请查看第 9.7 节）：
    • 调用 detslp_ll()，根据LLI检查接收机和基站观测数据是否有周跳；
    • 调用 detslp_gf()，利用几何无关组合进行周跳检测；
    • 调用 detslp_dop()，利用多普勒进行周跳检测（暂未使用）；
    • 更新 half-cycle, valid, flag。
  • 判断是否重置单差相位偏移（instantaneous 模式、中断次数超限或周跳）：
    • instantaneous 模式或中断次数超限：遍历每一颗卫星，若为 instantaneous 模式，或者卫星载波相位的中断次数 rtk->ssat[i-1].outc 大于配置中所设置的最大次数 rtk->opt.maxout，则重置载波偏移值为 0；
    • 遍历共视卫星，对单差相位偏移状态量的协方差阵加入过程噪声 rtk->opt.prn[0]*rtk->opt.prn[0]*fabs(tt)，如果发现有周跳或者异常值，则单差相位偏移状态量重置为0。
  • 更新相位偏移，遍历共视卫星：
    • 如果不是 IONOOPT_IFLC 消电离层组合：
      • 调用 sdobs() 函数，对基站、流动站观测值作差，求出载波相位单差观测值 cp、伪距单差观测值 pr；
      • 调用 sat2freq() 函数，获取载波频率；
      • 无效观测值，continue；
      • 相位偏差 bias[i] = 单差伪距 - (单差载波相位/光速)：bias[i]=cp-pr*freqi/CLIGHT;。
    • 如果是消电离层组合：获取两个频率的单差伪距，用消电离层组合的方法计算相位偏差 bias[i]；
    • 如果相位偏移不为0，有效卫星数 j++ ，offset 累加上与原来相位偏移的差。
  • 获取单差相位偏移平均值并更新 x：上一步累加有效卫星的相位偏差的差值到 offset，除以有效星的数 j，求得单差相位偏移平均值，所有有效卫星的相位偏差加上此平均值，完成相位偏移更新：rtk->x[IB(i,k,&rtk->opt)] += offset/j。
  • 单差相位偏移初始化：利用求得的 bias[i] 来对相位偏差参数为 0 的卫星，进行单差相位偏移状态量的初始化。

4. 源码注释

点击查看代码
c
static void udbias(rtk_t *rtk, double tt, const obsd_t *obs, const int *sat,
                   const int *iu, const int *ir, int ns, const nav_t *nav)
{
    double cp,pr,cp1,cp2,pr1,pr2,*bias,offset,freqi,freq1,freq2,C1,C2;
    int i,j,k,slip,reset,nf=NF(&rtk->opt);
    
    trace(3,"udbias  : tt=%.3f ns=%d\n",tt,ns);
    
    // 首先是对所有共视星进行周跳检测
    // 遍历每一颗卫星
    for (i=0;i<ns;i++) {
        // 调用 detslp_ll()，根据 LLI 检查接收机和基站观测数据是否有周跳
        /* detect cycle slip by LLI */  
        for (k=0;k<rtk->opt.nf;k++) rtk->ssat[sat[i]-1].slip[k]&=0xFC;
        detslp_ll(rtk,obs,iu[i],1);
        detslp_ll(rtk,obs,ir[i],2);

        // 调用 detslp_gf() 利用几何无关组合进行周跳检测
        /* detect cycle slip by geometry-free phase jump */ 
        detslp_gf(rtk,obs,iu[i],ir[i],nav);
        
        // 调用 detslp_dop()，利用多普勒进行周跳检测（但该函数由于时间跳变的原因，暂未使用）。
        /* detect cycle slip by doppler and phase difference */ //调用 detslp_dop 通过多普勒和相位差检查接收机和基站观测数据是否有周跳
        detslp_dop(rtk,obs,iu[i],1,nav);
        detslp_dop(rtk,obs,ir[i],2,nav);
        
        /* update half-cycle valid flag */
        for (k=0;k<nf;k++) {
            rtk->ssat[sat[i]-1].half[k]=
                !((obs[iu[i]].LLI[k]&2)||(obs[ir[i]].LLI[k]&2));
        }
    }

    // 遍历每一个频率
    for (k=0;k<nf;k++) {
        // 若为 instantaneous 模式，
        // 或者卫星载波相位的中断次数 rtk->ssat[i-1].outc 大于配置中所设置的最大次数 rtk->opt.maxout，
        // 重置载波偏移值
        /* reset phase-bias if instantaneous AR or expire obs outage counter */
        for (i=1;i<=MAXSAT;i++) {
            
            reset=++rtk->ssat[i-1].outc[k]>(uint32_t)rtk->opt.maxout;
            
            if (rtk->opt.modear==ARMODE_INST&&rtk->x[IB(i,k,&rtk->opt)]!=0.0) {
                initx(rtk,0.0,0.0,IB(i,k,&rtk->opt));
            }
            else if (reset&&rtk->x[IB(i,k,&rtk->opt)]!=0.0) {
                initx(rtk,0.0,0.0,IB(i,k,&rtk->opt));
                trace(3,"udbias : obs outage counter overflow (sat=%3d L%d n=%d)\n",
                      i,k+1,rtk->ssat[i-1].outc[k]);
                rtk->ssat[i-1].outc[k]=0;
            }
            if (rtk->opt.modear!=ARMODE_INST&&reset) {
                rtk->ssat[i-1].lock[k]=-rtk->opt.minlock;
            }
        }
        // 对共视星进行循环，对单差相位偏移状态量的协方差阵加入过程噪声 rtk->opt.prn[0]*rtk->opt.prn[0]*fabs(tt)
        // 如果发现有周跳或者异常值，则单差相位偏移状态量重置为0。
        /* reset phase-bias if detecting cycle slip */
        for (i=0;i<ns;i++) {
            j=IB(sat[i],k,&rtk->opt);
            rtk->P[j+j*rtk->nx]+=rtk->opt.prn[0]*rtk->opt.prn[0]*fabs(tt);
            slip=rtk->ssat[sat[i]-1].slip[k];
            if (rtk->opt.ionoopt==IONOOPT_IFLC) slip|=rtk->ssat[sat[i]-1].slip[1];
            if (rtk->opt.modear==ARMODE_INST||!(slip&1)) continue;
            rtk->x[j]=0.0;
            rtk->ssat[sat[i]-1].lock[k]=-rtk->opt.minlock;
        }
        bias=zeros(ns,1);


        // 利用单差伪距、单差载波相位计算一个 单差相位偏移 估计值，来对单差相位偏移状态量进行更新
        
        /* estimate approximate phase-bias by phase - code 
        对共视星进行循环，利用“单差伪距”和“单差载波相位”计算一个“单差相位偏移”估计值，
        来对单差相位偏移状态量进行更新。由于如果忽略伪距误差，那么伪距减去载波相位，
        则应该是载波相位（m）的相位偏移(m)，所以这里计算的是一个大致的相位偏移bias[i]。
        如果配置为无电离层组合，计算则按照无电离层组合的方式来计算。
        最后，仅计算所有有效星（单差相位偏移状态不为0）的 offset = sum of (bias - phase-bias)。
        其实就是计算每颗有效星 bias[i] 与单差相位偏移状态量的偏差，然后把所有有效星的这个偏差值加起来，
        之后会除以有效星的个数，最终就是求一个偏差平均值 rtk->com_bias。
        */
       //对共视星进行循环
        for (i=j=0,offset=0.0;i<ns;i++) {
            //不是消电离层模式
            if (rtk->opt.ionoopt!=IONOOPT_IFLC) {
                cp=sdobs(obs,iu[i],ir[i],k); /* cycle */       // 载波相位单差观测值
                pr=sdobs(obs,iu[i],ir[i],k+NFREQ);             // 伪距单差观测值
                freqi=sat2freq(sat[i],obs[iu[i]].code[k],nav); // 载波相位频率
                if (cp==0.0||pr==0.0||freqi==0.0) continue;
                
                bias[i]=cp-pr*freqi/CLIGHT; // 大致的相位偏移
            }
            // 消电离层模式
            else {
                cp1=sdobs(obs,iu[i],ir[i],0);
                cp2=sdobs(obs,iu[i],ir[i],1);
                pr1=sdobs(obs,iu[i],ir[i],NFREQ);
                pr2=sdobs(obs,iu[i],ir[i],NFREQ+1);
                freq1=sat2freq(sat[i],obs[iu[i]].code[0],nav);
                freq2=sat2freq(sat[i],obs[iu[i]].code[1],nav);
                if (cp1==0.0||cp2==0.0||pr1==0.0||pr2==0.0||freq1==0.0||freq2<=0.0) continue;
                
                C1= SQR(freq1)/(SQR(freq1)-SQR(freq2));
                C2=-SQR(freq2)/(SQR(freq1)-SQR(freq2));
                
                bias[i]=(C1*cp1*CLIGHT/freq1+C2*cp2*CLIGHT/freq2)-(C1*pr1+C2*pr2);
            }
            if (rtk->x[IB(sat[i],k,&rtk->opt)]!=0.0) {
                offset += bias[i]-rtk->x[IB(sat[i],k,&rtk->opt)]; //与原来状态量中单差模糊度的偏差的总和
                j++;
            }
        }
     
        /* correct phase-bias offset to enssure phase-code coherency */
        if (j>0) {
            for (i=1;i<=MAXSAT;i++) {
                if (rtk->x[IB(i,k,&rtk->opt)]!=0.0) rtk->x[IB(i,k,&rtk->opt)]+=offset/j;
            }
        }
        // 利用求得的 bias[i] 来对没有进行初始化的卫星，进行单差相位偏移状态量的初始化。
        /* set initial states of phase-bias */
        for (i=0;i<ns;i++) {
            if (bias[i]==0.0||rtk->x[IB(sat[i],k,&rtk->opt)]!=0.0) continue;
            initx(rtk,bias[i],SQR(rtk->opt.std[0]),IB(sat[i],k,&rtk->opt));
        }
        free(bias);
    }
}