本科生学习GNSS算法中级教程（二）- rtklib多系统多频单点定位算法 - tgd修正以及代码实现

最新推荐文章于 2024-06-22 00:26:05 发布

梧桐Fighting

最新推荐文章于 2024-06-22 00:26:05 发布

阅读量1.5k

点赞数 2

分类专栏： GNSS 文章标签：机器学习人工智能算法

本文链接：https://blog.csdn.net/dong20081991/article/details/126612068

版权

GNSS 专栏收录该内容

23 篇文章 84 订阅

订阅专栏

如何修正码偏差

上一节介绍了码偏差产生的原理，以及我们为何要修正它。

给了一篇参考论文，论文中有一节专门介绍如何对多系统多频点的码片差进行修正。下面给出论文中的截图：

从上图中的公式我们可以知道，gps的P1和P2的修正量都是Tgd乘以一个系数。而GLO和GAL系统和GPS类似一样的修正逻辑。

BDS系统由于是以B3作为卫星钟差的参考基准，所以其修正公式如下：

代码实现

rtklib中单点定位程序调用修正码偏差的函数为pntpos.c->prange()，为和原始区分，我们支持多频修正的函数取名为prange_mulfreq()，并多传入一个频率的index，用来指示计算哪个频率的修正后的伪距。

GPS系统的修正逻辑如下，注意查看我增加的注释：

if (sys == SYS_GPS || sys == SYS_QZS)
        {
            gamma = SQR(FREQ1 / FREQ2);
            b1 = gettgd(sat, nav, 0); /* TGD (m) */
            switch (obs->code[k])
            {
            /*for L1C, we need to calibrate the DCB between P1 and C1. but here we ignore it*/
            case CODE_L1C:
            case CODE_L1P:
                tgd = b1;
                break;
            /*for L2P, calibrate tgd according the paper*/
            case CODE_L2P:
            case CODE_L2W:
                tgd = b1 * gamma;
                break;
            /*no L5 tgd info in broadcast nav, give the default value 0*/
            case CODE_L5Q:
                tgd = 0.0;
                break;
            }
        }

BDS系统修正方法与其他有区别，将BDS系统的修正逻辑粘贴如下：

else if (sys == SYS_CMP)
        {   
            switch (obs->code[k])
            {
            case CODE_L2I:
                tgd = gettgd(sat, nav, 0);
                break;
            case CODE_L7I:
                tgd = gettgd(sat, nav, 1);
                break;
            /*for BDS sys, the reference frequency is 6I. that's why the correction in 6I is zeror*/
            case CODE_L6I:
                tgd = 0.0;
                break;
            case CODE_L1P:
                tgd = gettgd(sat, nav, 2);
                break;
            case CODE_L5P:
                tgd = gettgd(sat, nav, 3);
                break;
            }
        }

BDS2支持的信号为2I/7I/6I，一共三频，一般认为按顺序叫作北斗的第一、二、三频点。而BDS3则支持更多的频点，其播发了2I/6I，并且为了与GPS/GAL兼容互操作，播发了新的1C以及5P。所有的诸如2I/1P/5P的频点叫法，都来自于rinex的中的频率定义，详情请参考该文档。

最后贴上整个prange_mulfreq()函数。GLO修正逻辑部分未做完全考虑，可能存在错误。个人不太喜欢也不太建议使用GLO系统，现有的GPS+BDS+GAL卫星数已完全足够使用。

函数的前半部分其实保留了原函数的功能，即对无电离层组合观测值的修正。

整个函数

/* psendorange with code bias correction -------------------------------------*/
static double prange_mulfreq(const obsd_t *obs, const nav_t *nav, const prcopt_t *opt, const int k,
                             double *var)
{
    double P1, P2, gamma, b1, b2, P, tgd;
    int sat, sys;

    sat = obs->sat;
    sys = satsys(sat, NULL);
    P1 = obs->P[0];
    P2 = obs->P[1];
    *var = 0.0;

    if (P1 == 0.0 || (opt->ionoopt == IONOOPT_IFLC && P2 == 0.0))
        return 0.0;

    /* P1-C1,P2-C2 DCB correction */
    if (sys == SYS_GPS || sys == SYS_GLO)
    {
        if (obs->code[0] == CODE_L1C)
            P1 += nav->cbias[sat - 1][1]; /* C1->P1 */
        if (obs->code[1] == CODE_L2C)
            P2 += nav->cbias[sat - 1][2]; /* C2->P2 */
    }
    if (opt->ionoopt == IONOOPT_IFLC)
    { /* dual-frequency */

        if (sys == SYS_GPS || sys == SYS_QZS)
        { /* L1-L2,G1-G2 */
            gamma = SQR(FREQ1 / FREQ2);
            return (P2 - gamma * P1) / (1.0 - gamma);
        }
        else if (sys == SYS_GLO)
        { /* G1-G2 */
            gamma = SQR(FREQ1_GLO / FREQ2_GLO);
            return (P2 - gamma * P1) / (1.0 - gamma);
        }
        else if (sys == SYS_GAL)
        { /* E1-E5b */
            gamma = SQR(FREQ1 / FREQ7);
            if (getseleph(SYS_GAL))
            {                                                    /* F/NAV */
                P2 -= gettgd(sat, nav, 0) - gettgd(sat, nav, 1); /* BGD_E5aE5b */
            }
            return (P2 - gamma * P1) / (1.0 - gamma);
        }
        else if (sys == SYS_CMP)
        { /* B1-B2 */
            gamma = SQR(((obs->code[0] == CODE_L2I) ? FREQ1_CMP : FREQ1) / FREQ2_CMP);
            if (obs->code[0] == CODE_L2I)
                b1 = gettgd(sat, nav, 0); /* TGD_B1I */
            else if (obs->code[0] == CODE_L1P)
                b1 = gettgd(sat, nav, 2); /* TGD_B1Cp */
            else
                b1 = gettgd(sat, nav, 2) + gettgd(sat, nav, 4); /* TGD_B1Cp+ISC_B1Cd */
            b2 = gettgd(sat, nav, 1);                           /* TGD_B2I/B2bI (m) */
            return ((P2 - gamma * P1) - (b2 - gamma * b1)) / (1.0 - gamma);
        }
        else if (sys == SYS_IRN)
        { /* L5-S */
            gamma = SQR(FREQ5 / FREQ9);
            return (P2 - gamma * P1) / (1.0 - gamma);
        }
    }
    else
    { /* single-freq */
        if (k < 0 || k >= (NFREQ + NEXOBS))
        {
            return 0.0;
        }
        P = obs->P[k];
        *var = SQR(ERR_CBIAS);
        tgd = 0.0;
        if (sys == SYS_GPS || sys == SYS_QZS)
        {
            gamma = SQR(FREQ1 / FREQ2);
            b1 = gettgd(sat, nav, 0); /* TGD (m) */
            switch (obs->code[k])
            {
            /*for L1C, we need to calibrate the DCB between P1 and C1. but here we ignore it*/
            case CODE_L1C:
            case CODE_L1P:
                tgd = b1;
                break;
            /*for L2P, calibrate tgd according the paper*/
            case CODE_L2P:
            case CODE_L2W:
                tgd = b1 * gamma;
                break;
            /*no L5 tgd info in broadcast nav, give the default value 0*/
            case CODE_L5Q:
                tgd = 0.0;
                break;
            }
        }
        else if (sys == SYS_GLO)
        {   /*for GLO, the logic may has error*/
            gamma = SQR(FREQ1_GLO / FREQ2_GLO);
            b1 = gettgd(sat, nav, 0); /* -dtaun (m) */
            b1 = b1 / (gamma - 1.0);
            switch (obs->code[k])
            {
            case CODE_L1C:
            case CODE_L1P:
                tgd = b1;
                break;
            case CODE_L2C:
            case CODE_L2P:
                tgd = b1 * gamma;
                break;
            }
        }
        
        else if (sys == SYS_GAL)
        {   /*for GAL sys, there is two ephemerises*/
            if (getseleph(SYS_GAL))
                b1 = gettgd(sat, nav, 0); /* BGD_E1E5a */
            else
                b1 = gettgd(sat, nav, 1); /* BGD_E1E5b */
            switch (obs->code[k])
            {
            case CODE_L1C:
                tgd = b1;
                break;
            case CODE_L5Q:
                gamma = SQR(FREQ1 / FREQ5);
                tgd = b1 * gamma;
                break;
            case CODE_L7Q:
                gamma = SQR(FREQ1 / FREQ7);
                tgd = b1 * gamma;
                break;
            }
        }
        else if (sys == SYS_CMP)
        {   
            switch (obs->code[k])
            {
            case CODE_L2I:
                tgd = gettgd(sat, nav, 0);
                break;
            case CODE_L7I:
                tgd = gettgd(sat, nav, 1);
                break;
            /*for BDS sys, the reference frequency is 6I. that's why the correction in 6I is zeror*/
            case CODE_L6I:
                tgd = 0.0;
                break;
            case CODE_L1P:
                tgd = gettgd(sat, nav, 2);
                break;
            case CODE_L5P:
                tgd = gettgd(sat, nav, 3);
                break;
            }
        }
        return P - tgd;
    }
    return P1;
}

代码已上传到国内版git上（gitee，主要是网速快）。代码链接请在个人公众号回复 git 获取。在此也建议大家使用git工具clone代码，这样可以便于查看代码的更新历史，以及方便拉取最新的代码。

另外本节的代码在mulfreq-spp分支，注意分支切换。master分支是原始的rtklib代码。另外虽然mulfreq-spp分支有一些其他的改动，但暂时结果的正确性无法保证，我还在开发中。

另外一个，建议大家使用vscode来查看代码。git和vscode的使用请自行搜索博客或者bilibili。

下一节我会讲一下如何使用git查看代码的历史以及使用vscode查看每次提交的改动，这样方便大家理解代码的修改思路。