给 Cssrlib加注释（1）--有关卫星位置计算

本文链接：https://blog.csdn.net/m0_59951791/article/details/136384318

本文是一位导航工程专业学生分享其阅读GitHub开源项目CSSRLIB中PPP-RTK服务PPP-HAS部分的学习经历，关注PPP相关知识和改正数应用，详细解析了数据准备、文件加载、频率选择和数据处理过程，特别是卫星位置计算与改正数的使用和校正。

摘要由CSDN通过智能技术生成

写在前面：

本人是大三导航工程专业学生，想了解轨道钟差融合方面知识，决定从头开始阅读最新开源的cssrlib代码的PPP-HAS部分(GitHub - hirokawa/cssrlib: Python Toolkit for satellite-based open PPP/PPP-RTK services)，主要学习PPP相关知识以及改正数的应用，打算从头开始对数据处理部分进行阅读，争取做到代码逐行理解，由于本人学识有限，且第一次发文章，如有谬误还望大家指出，另外，有哪些没讲清楚的也希望大家指出，希望做出一份能让大家都读懂的文章！

那我们的故事开始

数据准备

代码的第一板块是处理文件的选择，ep设置为开始处理时刻，obs和nav文件换成需要处理的观测文件和导航文件，file_has为采用作者提供的解码程序(sbf_decode)解码得到的十六进制改正数文件,atx_file是天线文件

note：后两个板块是文件加载和频率选择，这里不进行探讨

数据处理

在采用HAS模式解算时，传进ppppos()里的是三个参数cs(HAS十六进制改正数，obs观测文件，nav星历文件)，接下来解算卫星位置(satposs)，obs表示精密轨道产品，在HAS里不需要。

satposs()函数--改正数与星历选择

初始化后进入大循环，首先提取卫星号，并判断系统：

sat = obs.sat[i]
sys, _ = sat2prn(sat)

在这个工程文件中，卫星的存储顺序为，GPS->GAL->QZS->BDS->GLO,与gamp相同的思路，根据卫星号判断系统(sat2prn)，假设卫星号为42，超出了GPS的32颗卫星(GALMIN)但是少于GAL和GPS的总卫星数68(QZSMIN)，那么卫星号变成10，sys识别为GAL(1)

下一步理解为为卫星选择伪距，由于解算顺序是按照观测文件的顺序来解算，因此伪距直接按顺序选择就好：

 pr = obs.P[i, 0]

这里[i,0]中的0表示选择第一频率的伪距，观测文件是双频，因此有0，1之分。

接下来是计算卫星发射时刻：

t = timeadd(obs.t, -pr/rCST.CLIGHT)

这里obs.t表示接收到卫星信号的时刻，并不是卫星信号发射时刻

nav.ephopt表示选择解算模式，4表示采用精密产品，这里进入else语句

紧接着判断是否有改正数：

cs.iodssr_c[sCType.ORBIT] == cs.iodssr

我对于这部分的理解是，这个工程里存储了12种类型的改正数(sCType)，iodssr_c是一个改正数有无的指示数组，如果有，则变成1，没有是-1

接下来找到卫星对应的下标：

idx = cs.sat_n.index(sat)

接下来提取版本号iode和轨道改正数dorb：

iode = cs.lc[0].iode[idx]
dorb = cs.lc[0].dorb[idx, :]

在判断id号没有改变后，说明轨道改正数和钟差改正数应该是在同一个位置(idx)，因此提取轨道改正数，如果有改变，会更新idx，方法同上

dclk = cs.lc[0].dclk[idx]

接下来判断改正数是否是数字（nan）

if np.isnan(dclk) or np.isnan(dorb@dorb):
    continue

接下来要选择星历，选择星历前注意B2b改正数用的是cnav，需要单独下载，不在广播星历中，而其余三系统使用的都是广播星历，作者在此处加了一个判断：

mode = cs.nav_mode[sys]

mode只有在sys=3，即北斗系统下是1，选择星历代码如下：

def findeph(nav, t, sat, iode=-1, mode=0):
    """ find ephemeric for sat """
    dt_p = 3600*4
    eph = None
    for eph_ in nav:
        if eph_.sat != sat:
            continue
        dt = timediff(t, eph_.toe)
        if (iode < 0 or eph_.iode == iode) and eph_.mode == mode and \
                abs(dt) < dt_p:
            eph = eph_
            break
    return eph

这段代码有以下几点注意事项：

这里通过遍历的方式寻找星历
计算接近时间dt时，用的是发射时间t，而不是接收时间
改正数的iode需要和星历的iode相匹配
只要小于阈值，就会进行赋值跳出循环，不会是最接近的时刻，会不会有问题？

接下来保存卫星健康信息和钟差

svh[i] = eph.svh
dt = eph2clk(t, eph)

eph2clk()函数通过广播星历计算卫星钟差，这里不做介绍，得到钟差后修正发射时间，让其更准确：

t = timeadd(t, -dt)

接下来计算卫星位置，由于我们的处理模式是HAS，本质上是先通过广播星历得到卫星位置，钟差，再根据改正数修正，所以这里直接通过广播星历计算位置（ECEF）、速度、钟差：

rs[i, :], vs[i, :], dts[i] = eph2pos(t, eph, True)

satposs()函数--卫星位置计算与改正数改正

通过eph2pos()函数计算广播星历下的卫星位置

 rs[i, :], vs[i, :], dts[i] = eph2pos(t, eph, True)

根据卫星系统选择地球常数

sys, prn = sat2prn(eph.sat)
    if sys == uGNSS.GAL:
        mu = rCST.MU_GAL
        omge = rCST.OMGE_GAL
    elif sys == uGNSS.BDS:
        mu = rCST.MU_BDS
        omge = rCST.OMGE_BDS
    else:  # GPS,QZS
        mu = rCST.MU_GPS
        omge = rCST.OMGE

进行时间的归一化[1]：

$\left\{\begin{array}{l} \mathrm{t}_k=\mathrm{t}_k-604800\left(\mathrm{t}_k>604800\right) \\ \mathrm{t}_k=\mathrm{t}_k+604800\left(\mathrm{t}_k \leq-604800\right) \end{array}\right.$

dt = dtadjust(t, eph.toe)
def dtadjust(t1, t2, tw=604800):
    """ calculate delta time considering week-rollover """
    dt = timediff(t1, t2)
    if dt > tw:
        dt -= tw
    elif dt < -tw:
        dt += tw
    return dt

再根据参考文献[1]中的公式计算卫星位置就好，代码是对应的，这里不一一列举

接下来通过改正数对卫星位置进行修正，通过接口文件计算坐标旋转矩阵[2]：

$\begin{gathered} e_t=\frac{\dot{x^s}}{\left|\dot{x}^s\right|} \\ e_w=\frac{x^s \times \dot{x}^s}{\left|x^s \times \dot{x}^s\right|} \\ \boldsymbol{e}_{\boldsymbol{n}}=\boldsymbol{e}_t \times \boldsymbol{e}_w \\ \boldsymbol{R}_{N T W}^{E C E F}=\left[\begin{array}{lll} \boldsymbol{e}_n & \boldsymbol{e}_t & \boldsymbol{e}_w \end{array}\right] \end{gathered}$

ea = vnorm(vs[i, :])
rc = np.cross(rs[i, :], vs[i, :])
ec = vnorm(rc)
er = np.cross(ea, ec)
A = np.array([er, ea, ec])

对改正数进行坐标系的旋转（从NTW到ECEF）

 dorb_e = dorb@A

接下来进行修正（符号与接口文件有差异）：

 rs[i, :] -= dorb_e
 dts[i] += dclk/rCST.CLIGHT

下面这段我也没咋看懂，但是结合作者的论文，推测是在计算SISRE

ers = vnorm(rs[i, :]-nav.x[0:3])
                dorb_ = -ers@dorb_e
                sis = dclk-dorb_
                if cs.lc[0].t0[1].time % 30 == 0 and \
                        timediff(cs.lc[0].t0[1], nav.time_p) > 0:
                    if abs(nav.sis[sat]) > 0:
                        nav.dsis[sat] = sis - nav.sis[sat]
                    nav.sis[sat] = sis

总结

改正数的选择直接采用的接收时刻进行索引
轨道改正数和钟差改正数的一致性通过mask_id控制
所谓改正数改正，直接在广播星历算好后加上改正数就行

qcedt()对于观测数据的预处理

在进行地球固体潮改正后，对所有卫星进行筛选，首先判断该颗卫星是否在观测文件中出现：

 if sat_i not in obs.sat:
            nav.edt[i, :] = 1
            continue

接下来找到观测文件中第一次出现该卫星的位置：

j = np.where(obs.sat == sat_i)[0][0]

找到这个的目的是为了判断该颗卫星有没有在前面算出轨道钟差信息，之前提到过，当时算卫星的顺序是根据观测文件来的，所以这里按照观测文件来给出索引值。

if np.isnan(rs[j, :]).any() or np.isnan(dts[j]):
    nav.edt[i, :] = 1
    if nav.monlevel > 0:
        nav.fout.write("{}  {} - edit - invalid eph\n"
                       .format(time2str(obs.t), sat2id(sat_i)))
    continue

接下来检查卫星健康信息：

if svh[j] > 0:
    nav.edt[i, :] = 1
    if nav.monlevel > 0:
        nav.fout.write("{}  {} - edit - satellite unhealthy\n"
                       .format(time2str(obs.t), sat2id(sat_i)))
    continue

检查卫星高度角，小于阈值则打上标记：

  _, e = gn.geodist(rs[j, :], rr_)
        _, el = gn.satazel(pos, e)
        if el < nav.elmin:
            nav.edt[i][:] = 1
            if nav.monlevel > 0:
                nav.fout.write("{}  {} - edit - low elevation {:5.1f} deg\n"
                               .format(time2str(obs.t), sat2id(sat_i),
                                       np.rad2deg(el)))
            continue

接下来通过观测文件判断是否失锁：

if obs.lli[j, f] == 1:
       nav.edt[i, f] = 1
       if nav.monlevel > 0:
             nav.fout.write("{}  {} - edit {:4s} - LLI\n"
             .format(time2str(obs.t), sat2id(sat_i),
               sigsCP[f].str()))
          continue

这里的lli指的是锁定丢失指示符（Loss of Lock Indicator）一个用来标记在接收到的卫星信号中相位连续性可能被中断的情况的指标，1表示有丢失，0表示没丢失，如果监测到有丢失，这个卫星的这个频率标记为1

接下来检查频率上的伪距和相位是不是0：

            if obs.P[j, f] == 0.0:
                nav.edt[i, f] = 1
                if nav.monlevel > 0:
                    nav.fout.write("{}  {} - edit {:4s} - invalid PR obs\n"
                                   .format(time2str(obs.t), sat2id(sat_i),
                                           sigsPR[f].str()))
                continue

            if obs.L[j, f] == 0.0:
                nav.edt[i, f] = 1
                if nav.monlevel > 0:
                    nav.fout.write("{}  {} - edit {:4s} - invalid CP obs\n"
                                   .format(time2str(obs.t), sat2id(sat_i),
                                           sigsCP[f].str()))
                continue

下面这段代码是根据信号类型来确定载噪比的阈值：

cnr_min = nav.cnr_min_gpy if sigsCN[f].isGPS_PY() else nav.cnr_min

对于 GPS P(Y) 信号，使用一个特定的阈值（nav.cnr_min_gpy）；而对于其他类型的信号，使用一个通用阈值（nav.cnr_min）

以上就做完了所有的预判断工作，主要的工作是给所有观测数据有问题的卫星及其频率打上标签,下一步进行滤波前的初始化。

参考文献

[1] 宋丹丹. 北斗导航定位解算算法的研究与软件实现[D].西安电子科技大学,2018.

[2] Galileo_HAS_SIS_ICD_v1.0