给 Cssrlib加注释（5）--星间单差非差非组合与矩阵构建代码与推导

最新推荐文章于 2024-08-04 22:34:46 发布

赵卷卷

最新推荐文章于 2024-08-04 22:34:46 发布

阅读量1.7k

点赞数 33

文章标签： python

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

版权

sdres()函数--单差处理与验前残差

首先根据系统进行循环，对于HAS服务就是GPS和Galileo，然后对于每个系统的不同频率循环两次，第一次针对伪距观测值，第二次针对相位观测值，在给卫星打上索引后，选取高度角最高的卫星为参考卫星：

  i = idx[np.argmax(el[idx])]

接下来对每颗卫星循环求解，首先获取电离层延迟的参考频率：

                    if sys == uGNSS.GLO:
                        freq0 = obs.sig[sys][uTYP.L][0].frequency(0)
                    else:
                        freq0 = obs.sig[sys][uTYP.L][0].frequency()

对于每个系统的每个信号类型都有自己的参考频率

下一步计算当前频率与参考频率的比值

                    if f < nf:  # carrier
                        sig = obs.sig[sys][uTYP.L][f]
                        if sys == uGNSS.GLO:
                            freq = sig.frequency(self.nav.glo_ch[sat[j]])
                        else:
                            freq = sig.frequency()
                        mu = -(freq0/freq)**2
                    else:  # code
                        sig = obs.sig[sys][uTYP.C][f % nf]
                        if sys == uGNSS.GLO:
                            freq = sig.frequency(self.nav.glo_ch[sat[j]])
                        else:
                            freq = sig.frequency()
                        mu = +(freq0/freq)**2

这里的if else语句判断是处理载波还是伪距，由于电离层延迟和频率的平方呈正比因此是平方项，对于载波相位观测值电离层延迟是负号，对于伪距观测值电离层延迟是正号

接下来跳过不可用的观测值：

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

跳过空值的观测：

                    if y[i, f] == 0.0 or y[j, f] == 0.0:
                        continue

跳过参考卫星：

                    if i == j:
                        continue

计算得到单差/双差观测值：

                    if mode == 0:  # DD
                        v[nv] = (y[i, f]-y[i+ns, f])-(y[j, f]-y[j+ns, f])
                    else:
                        #  Single-difference measurement
                        #
                        v[nv] = y[i, f] - y[j, f]

通过卫星到接收机的单位向量构建设计矩阵H：

H[nv, 0:3] = -e[i, :] + e[j, :]

这里设计矩阵的构建采用逐步扩展的方法，慢慢的和状态向量进行对应

对流层延迟采用参数估计的方式则会更新设计矩阵：

                    if self.nav.ntrop > 0:  # tropo is estimated

                        # SD troposphere
                        #
                        _, mapfwi = tropmapf(
                            obs.t, pos, el[i], model=self.nav.trpModel)
                        _, mapfwj = tropmapf(
                            obs.t, pos, el[j], model=self.nav.trpModel)

                        idx_i = self.IT(self.nav.na)
                        H[nv, idx_i] = mapfwi - mapfwj
                        v[nv] -= (mapfwi - mapfwj)*x[idx_i]

idx_x表示找到对流层延迟的位置，找到后更新设计矩阵，紧接着调整残差阵，x[idx_i]表示对流层的估计值

电离层延迟采用参数估计的方式也会更新设计矩阵：

                    if self.nav.niono > 0:  # iono is estimated
                        # SD ionosphere
                        #
                        idx_i = self.II(sat[i], self.nav.na)
                        idx_j = self.II(sat[j], self.nav.na)
                        H[nv, idx_i] = +mu
                        H[nv, idx_j] = -mu
                        v[nv] -= mu*(x[idx_i] - x[idx_j])

上述代码先找到了两颗卫星的在矩阵中的下标，并利用上文的比值更新设计矩阵，最后更新残差。

我觉得这里不做个详细的理论分析只看代码有点空洞，将会在下文对理论进行介绍并对争取公式进行详细推导。

接下来处理单差模糊度，同样的，先找到参考频率和当前频率所在的位置并计算波长：

                        idx_i = self.IB(sat[i], f, self.nav.na)
                        idx_j = self.IB(sat[j], f, self.nav.na)

                        if sys == uGNSS.GLO:
                            lami = sig.wavelength(self.nav.glo_ch[sat[i]])
                            lamj = sig.wavelength(self.nav.glo_ch[sat[j]])
                        else:
                            lami = sig.wavelength()
                            lamj = lami

接下来更新设计矩阵和残差：

                        H[nv, idx_i] = +lami
                        H[nv, idx_j] = -lamj
                        v[nv] -= lami*(x[idx_i] - x[idx_j])

接下来更新对于参考频率和当前频率的测量方差：

                        Ri[nv] = self.varerr(self.nav, el[i], f)
                        # measurement variance
                        Rj[nv] = self.varerr(self.nav, el[j], f)

这里应该是采用了高度角定权，最后标记了有效卫星：

                        self.nav.vsat[sat[i]-1, f] = 1
                        self.nav.vsat[sat[j]-1, f] = 1

对于伪距观测值则是根据高度角直接计算方差：

                    else:  # pseudorange

                        # measurement variance
                        Ri[nv] = self.varerr(self.nav, el[i], f)
                        # measurement variance
                        Rj[nv] = self.varerr(self.nav, el[j], f)

最后更新计数器即可

理论补充2--函数模型构建与公式推导（非差非组合）

可以看出cssrlib的函数模型是基于星间单差的非差非组合模型，周锋老师在他的博士论文里详细推导了非差非组合模型，因此在其基础上加入星间单差部分。

原始的GNSS观测模型可以表示为[6]：

$\begin{aligned} P_{r, j}^{s, Q} & =\rho_r^{s, Q}+c\left(d t_r-d t^{s, Q}\right)+\mathrm{Mw}_r^{s, Q} \cdot \mathrm{ZWD}_r+\gamma_j^Q \cdot I_{r, 1}^{s, Q} \\ & +c\left(d_{r, j}^{s, Q}-d_j^{s, Q}\right)+\varepsilon_{r, j}^{s, Q} \\ L_{r, j}^{s, Q} & \equiv \lambda_j^{s, Q} \cdot \Phi_{r, j}^{s, Q} \\ & =\rho_r^{s, Q}+c\left(d t_r-d t^{s, Q}\right)+\mathrm{Mw}_r^{s, Q} \cdot \mathrm{ZWD}_r-\gamma_j^Q \cdot I_{r, 1}^{s, Q} \\ & +\lambda_j^{s, Q}\left(N_{r, j}^{s, Q}+b_{r, j}^{s, Q}-b_j^{s, Q}\right)+\xi_{r, j}^{s, Q} \end{aligned}$

符号说明：

P：伪距

s：卫星号（PRN）

Q：卫星系统（G\R\E\J\C）

ρ：站星集合距离

dtr：接收机钟差

dts：卫星钟差

Mw：卫星高度角有关的湿投影函数

ZWD：测站天顶湿延迟

γ：频率相关的电离层延迟放大因子（仅与卫星系统有关，与卫星无关）

两个d：为频率相关的接收机和卫星端的非校正伪距硬件延迟（UCD）

两个b：为频率相关的接收机和卫星端的非校正相位硬件延迟（UPD）

I：是频率f1上的斜电离层延迟

N：是载波整周模糊度

还有两项是噪声、多路径误差与未模型化的误差

固体潮、相对论效应、对流层干延迟、天线相位缠绕等模型改正误差没写

周老师论文里引入了精密钟差产品，HAS服务也是基于双频观测值，因此和精密产品一样也会引入基准偏差，推导相同

定义

$\left\{\begin{array}{l} \alpha_{m n}^Q=\frac{\left(f_m^{s, Q}\right)^2}{\left(f_m^{s, Q}\right)^2-\left(f_n^{s, Q}\right)^2}, \beta_{m n}^Q=-\frac{\left(f_n^{s, Q}\right)^2}{\left(f_m^{s, Q}\right)^2-\left(f_n^{s, Q}\right)^2} \\ \mathrm{DCB}_{P_m P_n}^{s, Q}=d_m^{s, Q}-d_n^{s, Q}, \mathrm{DCB}_{r, P_m P_n}^{s, Q}=d_{r, m}^{s, Q}-d_{r, n}^{s, Q} \\ d_{\mathrm{IF} m n}^{s, Q}=\alpha_{m n}^Q \cdot d_m^{s, Q}+\beta_{m n}^Q \cdot d_n^{s, Q}, d_{r, \mathrm{~F}_{m n}}^{s, Q}=\alpha_{m n}^Q \cdot d_{r, m}^{s, Q}+\beta_{m n}^Q \cdot d_{r, n}^{s, Q} \end{array}\right.$

符号说明：

α、β：频率相关的放大因子

DCB：频率相关的卫星和接收机差分码偏差

第一个dIF：卫星UCD无电层组合第二个dIF：接收机UCD无电离层组合

其中，卫星端UCD和钟差线性相关，被钟差吸收，HAS服务提供的改正数也是基于无电离层组合：

$\begin{aligned} c d t_{\mathrm{IF}_{12}}^{s, Q} & =c\left[d t^{s, Q}+\left(\alpha_{12}^Q \cdot d_1^{s, Q}+\beta_{12}^Q \cdot d_2^{s, Q}\right)+d D^Q\right] \\ & =c\left(d t^{s, Q}+d_{\mathrm{IF}_{12}}^{s, Q}+d D^Q\right) \end{aligned}$

dD是卫星钟差估计引入的参考基准引起的偏差，被接收机钟差吸收。

把钟差改正数带入原始式子并进行线性化：

$\begin{aligned} p_{r, j}^{s, Q}= & \mathbf{u}_r^{s, Q} \cdot \mathbf{x}+c d t_r^Q+c d_{r, j}^{s, Q}+\mathrm{Mw}_r^{s, Q} \cdot \mathrm{ZWD}_r +\gamma_j^Q \cdot\left(I_{r, 1}^{s, Q}-\beta_{12}^Q \cdot c \mathrm{DCB}_{P_1 P_2}^{s, Q}\right)+\varepsilon_{r, j}^{s, Q} \\ l_{r, j}^{s, Q} & =\mathbf{u}_r^{s, Q} \cdot \mathbf{x}+c d t_r^Q+c d_{\mathrm{IF}_{12}}^{s, Q}+\mathrm{Mw}_r^{s, Q} \cdot \mathrm{ZWD}_r -\gamma_j^Q \cdot I_{r, 1}^{s, Q}+\lambda_j^{s, Q}\left(N_{r, j}^{s, Q}+b_{r, j}^{s, Q}-b_j^{s, Q}\right)+\xi_{r, j}^{s, Q} \end{aligned}$

符号说明：

p,l：表示伪距和载波的观测量减去计算量（就是原始观测值减去展开后的初值）

u：表示方向余弦

x：表示三维增量

接收机端伪距硬件延迟可同时被接收机钟差和电离层参数吸收：

$\left\{\begin{array} { l } { c d _ { r , 1 } ^ { Q } = a + \gamma _ { 1 } ^ { Q } \cdot b } \\ { c d _ { r , 2 } ^ { Q } = a + \gamma _ { 2 } ^ { Q } \cdot b } \end{array} \Rightarrow \left\{\begin{array}{l} a=\alpha_{12}^Q \cdot c d_{r, 1}^Q+\beta_{12}^Q \cdot c d_{r, 2}^Q=c d_{r, \mathrm{IF}_{12}}^Q \\ b=\frac{1}{1-\gamma_2^Q} \cdot\left(c d_{r, 1}^Q-c d_{r, 2}^Q\right)=\beta_{12}^Q \cdot c \mathrm{DCB}_{r, P_1 P_2}^Q \end{array}\right.\right.$

重新整理两式可得：

$\begin{aligned} p_{r, j}^{s, Q}= & \mathbf{u}_r^{s, Q} \cdot \mathbf{x}+c d \bar{t}_r^Q+\mathrm{Mw}_r^{s, Q} \cdot \mathrm{ZWD}_r \\ & +\gamma_j^Q \cdot \bar{I}_{r, 1}^{s, Q}+\varepsilon_{r, j}^{s, Q} \\ l_{r, j}^{s, Q} & =\mathbf{u}_r^{s, Q} \cdot \mathbf{x}+c d \bar{t}_r^Q+\mathrm{Mw}_r^{s, Q} \cdot \mathrm{ZWD}_r \\ & -\gamma_j^Q \cdot \bar{I}_{r, 1}^{s, Q}+\bar{N}_{r, j}^{s, Q}+\xi_{r, j}^{s, Q} \end{aligned}$

$\left\{\begin{aligned} & c d \bar{t}_r^Q=c d t_r^Q+c d_{r, \mathrm{FF}_{12}}^Q \\ & \bar{I}_{r, 1}^{s, Q}= I_{r, 1}^{s, Q}+c \beta_{12}^Q\left(\mathrm{DCB}_{r, P_1 P_2}^Q-\mathrm{DCB}_{P_1 P_2}^{s, Q}\right) \\ & \bar{N}_{r, j}^{s, Q}= \lambda_j^{s, Q}\left(N_{r, j}^{s, Q}+b_{r, j}^{s, Q}-b_j^{s, Q}\right)+c\left(d_{\mathrm{IF}_{12}}^{s, Q}-d_{r, \mathrm{IF}_{12}}^Q\right) \\ & \quad+\frac{c \gamma_j^Q}{1-\gamma_2^Q}\left(\mathrm{DCB}_{r, P_1 P_2}^Q-\mathrm{DCB}_{P_1 P_2}^{s, Q}\right) \end{aligned}\right.$

HAS服务目前只应用于GPS和Galileo系统，这里不考虑格洛纳斯系统，以上的公式都是周老师博士论文里的公式[6]，下面推导星间单差：

首先是伪距：

接下来载波：

其实就是少了接收机钟差项，接下来可知，需要估计的参数有：三维坐标增量x，对流层湿延迟ZWD，电离层延迟I，模糊度N（两个频率的）总数为4+3n，因此设计矩阵的构建如下：

设计矩阵的大小为4m*n

参考文献

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

[2] Galileo_HAS_SIS_ICD_v1.0

[3] 李昕. 多频率多星座GNSS快速精密定位关键技术研究[D].武汉大学,2022.

[4] 《GNSS精密单点定位理论方法及其应用》

[5] 李星星.GNSS精密单点定位及非差模糊度快速确定方法研究[D].武汉大学,2013.

[6] 周锋.多系统GNSS非差非组合精密单点定位相关理论和方法研究[D].华东师范大学,2018.

赵卷卷

关注

33
点赞
踩
24

收藏

觉得还不错? 一键收藏
1
评论
给 Cssrlib加注释（5）--星间单差非差非组合与矩阵构建代码与推导

这里的if else语句判断是处理载波还是伪距，由于电离层延迟和频率的平方呈正比因此是平方项，对于载波相位观测值电离层延迟是负号，对于伪距观测值电离层延迟是正号。可以看出cssrlib的函数模型是基于星间单差的非差非组合模型，周锋老师在他的博士论文里详细推导了非差非组合模型，因此在其基础上加入星间单差部分。idx_x表示找到对流层延迟的位置，找到后更新设计矩阵，紧接着调整残差阵，x[idx_i]表示对流层的估计值。上述代码先找到了两颗卫星的在矩阵中的下标，并利用上文的比值更新设计矩阵，最后更新残差。
复制链接

扫一扫