RTKLIB算法优化之北斗伪距多路径代码

不同于GPS，BDS伪距观测值中存在与高度角强相关的系统性偏差。该现象最早由Hauschild 等人发现(Hauschild et al. 2012)，随后Gisbert 等(Gisbert et al. 2012)，Montenbruck 等(Montenbruck et al.2013b)都对北斗伪距偏差做了进一步的研究，表明北斗伪距偏差是与高度角及频率相关，而与接收机类型和方位角无关的系统偏差，并推断这种误差可能来自于BDS卫星端伪距多路径效应。对于BDS IGSO和MEO卫星该偏差变化幅度可达1m,严重影响了利用伪距观测信息的标准单点定位、单频PPP和双频PPP解，同样也势必会严重影响基于MW组合解算的宽巷模糊度的精度与稳定性,降低 BDSWLFCB估计的质量。2015年，Wanninger和 Beer利用大量的实测数据，建立了北斗卫星伪距偏差改正的经验公式，并验证了进行伪距偏差改正对单频精密单点定位精度的提高有重要的作用(Wanninger and Beer 2015)。

 BDmultipathCorr(&rtk, obs, n);

/* mutipath correct-------------------------------------------------------------
* BeiDou satellite-induced code pseudorange variations correct
* args   : rtk_t *rtk       IO  rtk control/result struct
           obsd_t *obs      IO  observation data
           int    n         I   number of observation data
           nav_t  *nav      I   navigation messages
* note   :
*
* -----------------------------------------------------------------------------*/
extern void BDmultipathCorr(rtk_t* rtk, obsd_t* obs, int n)
{
    int i, j, sat, b,prn,sys;
    double dmp[3], elev, a;
    const static double IGSOCOEF[3][10] = {        /* m */
        {-0.55,-0.40,-0.34,-0.23,-0.15,-0.04,0.09,0.19,0.27,0.35},    //B1
        {-0.71,-0.36,-0.33,-0.19,-0.14,-0.03,0.08,0.17,0.24,0.33},    //B2
        {-0.27,-0.23,-0.21,-0.15,-0.11,-0.04,0.05,0.14,0.19,0.32},    //B3
    };
    const static double MEOCOEF[3][10] = {        /* m */
        {-0.47,-0.38,-0.32,-0.23,-0.11,0.06,0.34,0.69,0.97,1.05},    //B1
        {-0.40,-0.31,-0.26,-0.18,-0.06,0.09,0.28,0.48,0.64,0.69},    //B2
        {-0.22,-0.15,-0.13,-0.10,-0.04,0.05,0.14,0.27,0.36,0.47},    //B3
    };

    for (i = 0; i < n && i < MAXOBS; i++) {
        sat = obs[i].sat;
        sys=satsys(sat,&prn);
        if (sys != SYS_CMP) continue;

        if (prn <= 5) continue;

        elev = rtk->ssat[sat - 1].azel[1] * R2D;

        if (elev <= 0.0) continue;

        for (j = 0; j < 3; j++) dmp[j] = 0.0;

        a = elev * 0.1;
        b = (int)a;

        if (prn >= 6 && prn < 11) { // IGSO(C06, C07, C08, C09, C10)
            if (b < 0) {
                for (j = 0; j < 3; j++) dmp[j] = IGSOCOEF[j][0];
            }
            else if (b >= 9) {
                for (j = 0; j < 3; j++) dmp[j] = IGSOCOEF[j][9];
            }
            else {
                for (j = 0; j < 3; j++) dmp[j] = IGSOCOEF[j][b] * (1.0 - a + b) + IGSOCOEF[j][b + 1] * (a - b);
            }
        }
        else if (prn >= 11) {   // MEO(C11, C12, C13, C14)
            if (b < 0) {
                for (j = 0; j < 3; j++) dmp[j] = MEOCOEF[j][0];
            }
            else if (b >= 9) {
                for (j = 0; j < 3; j++) dmp[j] = MEOCOEF[j][9];
            }
            else {
                for (j = 0; j < 3; j++) dmp[j] = MEOCOEF[j][b] * (1.0 - a + b) + MEOCOEF[j][b + 1] * (a - b);
            }
        }

        for (j = 0; j < 3; j++) obs[i].P[j] += dmp[j];
    }
}