给 Cssrlib加注释（4）--滤波前的误差项改正

数据处理

zdres()函数--非差分项误差改正

首先经过地球固体潮改正：

 if self.nav.tidecorr == uTideModel.SIMPLE:
            pos = ecef2pos(rr_)
            disp = tidedisp(gpst2utc(obs.t), pos)
        elif self.nav.tidecorr == uTideModel.IERS2010:
            pos = ecef2pos(rr_)
            disp = tidedispIERS2010(gpst2utc(obs.t), pos)
        else:
            disp = np.zeros(3)
        rr_ += disp

然后计算对流层干延迟和对流层湿延迟：

        trop_hs, trop_wet, _ = tropmodel(obs.t, pos,
                                         model=self.nav.trpModel)

如果有对流层改正数或者电离层改正数的话：

 if self.nav.trop_opt == 2 or self.nav.iono_opt == 2:  # from cssr
            inet = cs.find_grid_index(pos)
            dlat, dlon = cs.get_dpos(pos)
        else:
            inet = -1
​
        if self.nav.trop_opt == 2:  # trop from cssrz
            trph, trpw = cs.get_trop(dlat, dlon)
            trop_hs0, trop_wet0, _ = tropmodel(obs.t, [pos[0], pos[1], 0],
                                               model=self.nav.trpModel)
            r_hs = trop_hs/trop_hs0
            r_wet = trop_wet/trop_wet0
​
        if self.nav.iono_opt == 2:  # iono from cssr
            stec = cs.get_stec(dlat, dlon)

由于HAS服务不提供这两种改正数，因此这里不进行探讨。

首先跳过无效的卫星：

 if np.any(self.nav.edt[sat-1, :] > 0):
                continue

其次跳过没有改正数的：

            if inet > 0 and sat not in cs.lc[inet].sat_n:
                continue

inet表示找不到改正数，第二个判断应该是改正数信息没有包含当前卫星。

提取信号，波长，与频率：

            sigsPR = obs.sig[sys][uTYP.C]
            sigsCP = obs.sig[sys][uTYP.L]
​
            # Wavelength
            #
            if sys == uGNSS.GLO:
                lam = np.array([s.wavelength(self.nav.glo_ch[sat])
                                for s in sigsCP])
                frq = np.array([s.frequency(self.nav.glo_ch[sat])
                               for s in sigsCP])
            else:
                lam = np.array([s.wavelength() for s in sigsCP])
                frq = np.array([s.frequency() for s in sigsCP])

紧接着进入相位偏差与伪距偏差改正(DCB)，由于只采用广播星历，这里进入else语句，首先进行相位偏差修正：

                if cs.lc[0].cstat & (1 << sCType.CBIAS) == (1 << sCType.CBIAS):
                    cbias = self.find_bias(cs, sigsPR, sat)
​
                if inet > 0 and cs.lc[inet].cstat & (1 << sCType.CBIAS) == \
                        (1 << sCType.CBIAS):
                    cbias += self.find_bias(cs, sigsPR, sat, inet)

我们逐行分析下这段代码：

首先cstat是一个状态标志位，表示当前时刻改正数信息是否有效，cstate & 1 << sCType.CBIAS表示，如果cstate等于1（改正数有效），那么应该还是1 << sCType.CBIAS，如果相等就找改正数，找相位改正数的代码如下：

首先进行初始化和检查是否有信号：

    def find_bias(self, cs, sigref, sat, inet=0):
        """ find satellite signal bias from correction """
        nf = len(sigref)
        v = np.zeros(nf)
​
        if nf == 0:
            return v

提取第一频率的改正数：

        ctype = sigref[0].typ
        if ctype == uTYP.C:
            if cs.lc[inet].cbias is None or \
                    sat not in cs.lc[inet].cbias.keys():
                return v
            sigc = cs.lc[inet].cbias[sat]
        else:
            if cs.lc[inet].pbias is None or \
                    sat not in cs.lc[inet].pbias.keys():
                return v
            sigc = cs.lc[inet].pbias[sat]

if通过代码提取相位偏差，不通过提取伪距偏差，在有改正数且改正数包含该卫星的情况下。

如果是HAS服务，还会进一步判断：

        if cs.cssrmode in [sc.GAL_HAS_SIS, sc.GAL_HAS_IDD]:
            if ctype == uTYP.C and rSigRnx('GC2P') in sigc.keys():
                sigc[rSigRnx('GC2W')] = sigc[rSigRnx('GC2P')]
            if ctype == uTYP.L and rSigRnx('GL2P') in sigc.keys():
                sigc[rSigRnx('GL2W')] = sigc[rSigRnx('GL2P')]

这里的意思是即如果信号类型为C或L，且存在特定的信号（如GC2P），则会将这个信号的偏置复制给另一个信号（如GC2W），根据接口文件[2]的表20可以看到，并没有提供GPS系统C2W和L2W的改正数，推测是这两个信道的改正数和C2P,L2P相同。

最后，函数遍历所有信号参考，如果信号参考在偏置数据sigc中，它会直接将对应的偏置值赋给向量v，如果直接的信号没有找到，它会尝试查找转换为另一属性（如X）的信号，并将对应的偏置值赋给v：

        for k, sig in enumerate(sigref):
            if sig in sigc.keys():
                v[k] = sigc[sig]
            elif sig.toAtt('X') in sigc.keys():
                v[k] = sigc[sig.toAtt('X')]

至于为什么会找两遍，我现在也不是很清楚，有知道的小伙伴可以在评论区讨论下。

相同的伪距偏差也是相同的提取方法，接下来判断改正数类型，符号是否相反需要参考对应的ICD接口文件：

                if cs.cssrmode in [sc.QZS_CLAS, sc.BDS_PPP, sc.PVS_PPP]:
                    pbias = -pbias
                    cbias = -cbias

在排除无效的相位改正数后，接下来进行地球自转改正并计算卫星高度角：

            r, e[i, :] = geodist(rs[i, :], rr_)
            _, el[i] = satazel(pos, e[i, :])
            if el[i] < self.nav.elmin:
                continue

其中，e表示接收机与卫星的单位矢量、el表示卫星高度角

接下来进行相对论效应的改正：

 relatv = shapiro(rs[i, :], rr_)

接下来计算对流层映射函数：

           mapfh, mapfw = tropmapf(obs.t, pos, el[i],
                     model=self.nav.trpModel)

这部分的计算是用于对流层改正数的使用：

            # Tropospheric delay
            #
            if self.nav.iono_opt == 2:  # from cssr
                trop = mapfh*trph*r_hs+mapfw*trpw*r_wet
            else:
                trop = mapfh*trop_hs + mapfw*trop_wet

紧接着计算电离层延迟：

            # Ionospheric delay
            if self.nav.iono_opt == 2:  # from cssr
                idx_l = cs.lc[inet].sat_n.index(sat)
                iono = np.array([40.3e16/(f*f)*stec[idx_l] for f in frq])
            else:
                iono = np.zeros(nf)

这里没有改正数初始化为零矩阵说明在后面可能采用参数估计的方式解算电离层延迟。

接下来进行天线相位缠绕改正：

            if self.nav.phw_opt > 0:
                phw_mode = (False if self.nav.phw_opt == 2 else True)
                self.nav.phw[sat-1] = windupcorr(obs.t, rs[i, :], vs[i, :],
                                                 rr_, self.nav.phw[sat-1],
                                                 full=phw_mode)
​
                # cycle -> m
                phw = lam*self.nav.phw[sat-1]
            else:
                phw = np.zeros(nf)

接下来是频率标注，

            sig0 = None
            if cs is not None:
​
                if cs.cssrmode in (sc.GAL_HAS_SIS, sc.GAL_HAS_IDD,
                                   sc.QZS_MADOCA):
​
                    if sys == uGNSS.GPS:
                        sig0 = (rSigRnx("GC1W"), rSigRnx("GC2W"))
                    elif sys == uGNSS.GAL:
                        sig0 = (rSigRnx("EC1C"), rSigRnx("EC7Q"))
                    elif sys == uGNSS.QZS:
                        sig0 = (rSigRnx("JC1C"), rSigRnx("JC2S"))
                    elif sys == uGNSS.GLO:
                        sig0 = (rSigRnx("RC1C"), rSigRnx("RC2C"))
​
                elif cs.cssrmode == sc.BDS_PPP:
​
                    if sys == uGNSS.GPS:
                        sig0 = (rSigRnx("GC1W"), rSigRnx("GC2W"))
                    elif sys == uGNSS.BDS:
                        sig0 = (rSigRnx("CC6I"),)

不同的改正数是基于不同的频率，比如HAS服务是基于C1W,C2W，C1C和C7Q频率的，具体改正数的对应频率需要参考ICD。

接下来进行天线相位中心改正：

            # Receiver/satellite antenna offset
            #
            antrPR = antModelRx(self.nav, pos, e[i, :], sigsPR, rtype)
            antrCP = antModelRx(self.nav, pos, e[i, :], sigsCP, rtype)
​
            if self.nav.ephopt == 4:
​
                antsPR = antModelTx(
                    self.nav, e[i, :], sigsPR, sat, obs.t, rs[i, :])
                antsCP = antModelTx(
                    self.nav, e[i, :], sigsCP, sat, obs.t, rs[i, :])
​
            elif cs is not None and (cs.cssrmode == sc.GAL_HAS_SIS or
                                     cs.cssrmode == sc.GAL_HAS_IDD or
                                     cs.cssrmode == sc.QZS_MADOCA or
                                     cs.cssrmode == sc.BDS_PPP):
​
                antsPR = antModelTx(self.nav, e[i, :], sigsPR,
                                    sat, obs.t, rs[i, :], sig0)
                antsCP = antModelTx(self.nav, e[i, :], sigsCP,
                                    sat, obs.t, rs[i, :], sig0)
​
            else:
​
                antsPR = [0.0 for _ in sigsPR]
                antsCP = [0.0 for _ in sigsCP]

注意对于接收机端和卫星端是不一样的，频率的选择方面有所区别，在没有cssr改正数时频率为none

接下来是对上一时刻计算得到的误差改正的总结与归纳：

            prc[i, :] = trop + antrPR + antsPR + iono - cbias
            cpc[i, :] = trop + antrCP + antsCP - iono - pbias + phw
            r += relatv - _c*dts[i]

接下来按照频率，进行观测向量y的构建：

            for f in range(nf):
                y[i, f] = obs.L[i, f]*lam[f]-(r+cpc[i, f])
                y[i, f+nf] = obs.P[i, f]-(r+prc[i, f])

y的维数应该是：i * (nf*2) i表示观测到的卫星数量，nf表示频率数量。