目前，多星座GNSS并存与发展的局面已初步形成，多模GNSS兼容与互操作以及多频多模GNSS数据融合已成为GNSS导航定位领域的主要发展趋势^[1-2]。相对于单系统，多系统可以显著增加可见卫星个数，改善卫星空间几何结构，提高系统抵御粗差的能力以及完备性。

如何合理地确定多系统观测值的权比是使用多卫星导航系统的关键之一。由于多系统组合定位中各系统之间的噪声不同，轨道计算精度略有差异，如果各类观测值之间权比不合理，会导致验后单位权方差有偏，影响参数估值的最优性^[3]。利用先验精度确定的观测值权矩阵通常不能完全反映各类观测值的实际精度或精度结果可信度不够。为了得到较好的组合定位结果，应尽可能合理地确定观测值的权。Helmert方差分量估计是通过迭代计算自适应确定不同类观测值权比的验后方差分量估计方法，在数据处理领域得到广泛应用，并取得了较好的效果^[4-7]。但是，对于GPS/GLONASS/BDS多系统组合定位中观测值权比确定尚缺乏深入研究。基于此，本文给出了GPS/GLONASS/BDS组合单点定位和基线解算模型，并引入Helmert方差分量估计，对不同系统观测值合理定权。实验结果表明，在GPS/GLONASS/BDS组合定位中，经过Helmert方差分量估计对不同系统观测值进行合理定权，可以得到各系统各类观测值的最佳权比，提高GPS/GLONASS/BDS组合定位的精度。

1 多星座GNSS组合单点定位模型

单系统单点定位误差方程为：

$ \mathit{\boldsymbol{V}}_p^{\rm{*}} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{B}}^{\rm{*}}}}&{-{\mathit{\boldsymbol{e}}^{\rm{*}}}} \end{array}} \right]\left[\begin{array}{c} \mathit{\boldsymbol{X}}\\ d{T^{\rm{*}}} \end{array} \right] - \mathit{\boldsymbol{L}}_p^{\rm{*}} $

(1)

式中，V为观测值残差向量，上标*=G、R、C，分别为GPS、GLONASS、BDS系统(下同)，下标p为伪距观测值，B为设计矩阵，元素为接收机与卫星视线方向单差单位向量，e为元素均为1的列矩阵，X为坐标参数，dT为接收机钟差，L为观测向量。

根据卫星高度角对GPS、GLONASS、BDS伪距观测值分别进行初始定权：

$ {\mathit{\boldsymbol{P}}_p}(i, i) = {(\sin ({E_i}))^2} $

(2)

对应的法方程为：

$ \left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{N}}_{11, p}^{\rm{*}}}&{\mathit{\boldsymbol{N}}_{12, p}^{\rm{*}}}\\ {\mathit{\boldsymbol{N}}_{21, p}^{\rm{*}}}&{\mathit{\boldsymbol{N}}_{22, p}^{\rm{*}}} \end{array}} \right]\left[\begin{array}{c} {\mathit{\boldsymbol{\hat X}}}\\ d{{\hat T}^{\rm{*}}} \end{array} \right] = \left[\begin{array}{c} \mathit{\boldsymbol{W}}_{1, p}^{\rm{*}}\\ \mathit{\boldsymbol{W}}_{2, p}^{\rm{*}} \end{array} \right] $

(3)

消除接收机钟差参数后，法方程变为：

$ $\begin{array}{l} (\mathit{\boldsymbol{N}}_{11, p}^{\rm{*}} - \mathit{\boldsymbol{N}}_{12, p}^{\rm{*}}{(\mathit{\boldsymbol{N}}_{22, p}^{\rm{*}})^{ - 1}}\mathit{\boldsymbol{N}}_{21, p}^{\rm{*}})\mathit{\boldsymbol{\hat X}} - \\ ({\mathit{\boldsymbol{W}}}_{1, p}^{\rm{*}} - {\mathit{\boldsymbol{N}}}_{12, p}^{\rm{*}}{({\mathit{\boldsymbol{N}}}_{22, p}^{\rm{*}})^{ - 1}}{\mathit{\boldsymbol{W}}}_{2, p}^{\rm{*}})\\ = {\mathit{\boldsymbol{N}}}_p^{\rm{*}}{\mathit{\boldsymbol{\hat X}}} - {\mathit{\boldsymbol{W}}}_p^{\rm{*}} = 0 \end{array}$ $

(4)

叠加法方程，求得测站三维坐标数为：

$ \begin{array}{l} {\mathit{\boldsymbol{\hat X}}} = {({\mathit{\boldsymbol{N}}}_p^{\rm{G}} + {\mathit{\boldsymbol{N}}}_p^{\rm{R}} + {\mathit{\boldsymbol{N}}}_p^{\rm{C}})^{-1}}({\mathit{\boldsymbol{W}}}_p^{\rm{G}} + {\mathit{\boldsymbol{W}}}_p^{\rm{R}} + {\mathit{\boldsymbol{W}}}_p^{\rm{C}})\\ = {({{\mathit{\boldsymbol{N}}}_p})^{-1}}{{\mathit{\boldsymbol{W}}}_p} \end{array} $

(5)

2 多星座GNSS组合基线解算模型

单系统基线解算误差方程为：

$ {\mathit{\boldsymbol{V}}}_\varphi ^{\rm{*}} = \left[{\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{B}}}^{\rm{*}}}}&{{{\mathit{\boldsymbol{C}}}^{\rm{*}}}} \end{array}} \right]\left[\begin{array}{c} {\mathit{\boldsymbol{X}}}\\ {{\mathit{\boldsymbol{n}}}^{\rm{*}}} \end{array} \right] -{\mathit{\boldsymbol{L}}}_\varphi ^{\rm{*}} $

(6)

式中，下标φ为载波相位观测值，C为设计矩阵，n为双差模糊度参数。

根据高度角模型对GPS、GLONASS、BDS载波相位观测值分别进行初始定权:

$ {\mathit{\boldsymbol{P}}}(i, i) = \frac{1}{{1 + \frac{{\sin {E_{{\rm{ref}}}}}}{{\sin {E_i}}}}} $

(7)

式中，E_ref、E_i分别为参考卫星、非参考卫星的高度角。

对应的法方程为：

$ \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{N}}}_{11, \varphi }^{\rm{*}}}&{{\mathit{\boldsymbol{N}}}_{12, \varphi }^{\rm{*}}}\\ {{\mathit{\boldsymbol{N}}}_{21, \varphi }^{\rm{*}}}&{{\mathit{\boldsymbol{N}}}_{22, \varphi }^{\rm{*}}} \end{array}} \right]\left[\begin{array}{c} {{\mathit{\boldsymbol{\hat X}}}}\\ {{{\mathit{\boldsymbol{\hat n}}}}^{\rm{*}}} \end{array} \right] = \left[\begin{array}{c} {\mathit{\boldsymbol{W}}}_{1, \varphi }^{\rm{*}}\\ {\mathit{\boldsymbol{W}}}_{2, \varphi }^{\rm{*}} \end{array} \right] $

(8)

消除模糊度参数后，法方程变为：

$ \begin{array}{l} ({\mathit{\boldsymbol{N}}}_{11, \varphi }^{\rm{*}} - {\mathit{\boldsymbol{N}}}_{12, \varphi }^{\rm{*}}{({\mathit{\boldsymbol{N}}}_{22, \varphi }^{\rm{*}})^{ - 1}}{\mathit{\boldsymbol{N}}}_{21, \varphi }^{\rm{*}}){\mathit{\boldsymbol{\hat X}}} - ({\mathit{\boldsymbol{W}}}_{1, \varphi }^{\rm{*}} - {\mathit{\boldsymbol{N}}}_{12, \varphi }^{\rm{*}}{({\mathit{\boldsymbol{N}}}_{22, \varphi }^{\rm{*}})^{ - 1}}{\mathit{\boldsymbol{W}}}_{2, \varphi }^{\rm{*}})\\ = {\mathit{\boldsymbol{N}}}_\varphi ^{\rm{*}}{\mathit{\boldsymbol{\hat X}}} - {\mathit{\boldsymbol{W}}}_\varphi ^{\rm{*}} = 0 \end{array} $

(9)

叠加法方程，求得测站三维坐标改正数为：

$ $\begin{array}{l} {\mathit{\boldsymbol{\hat X}}} = {({\mathit{\boldsymbol{N}}}_\varphi ^{\rm{G}} + {\mathit{\boldsymbol{N}}}_\varphi ^{\rm{R}} + {\mathit{\boldsymbol{N}}}_\varphi ^{\rm{C}})^{ - 1}}({\mathit{\boldsymbol{W}}}_\varphi ^{\rm{G}} + {\mathit{\boldsymbol{W}}}_\varphi ^{\rm{R}} + {\mathit{\boldsymbol{W}}}_\varphi ^{\rm{C}})\\ = {({{\mathit{\boldsymbol{N}}}_\varphi })^{ - 1}}{{\mathit{\boldsymbol{W}}}_\varphi } \end{array}$ $

(10)

3 Helmert方差分量估计

GPS/GLONASS/BDS的Helmert方差分量估计过程为^[8]：

1) 进行最小二乘平差，求得(V^*)^TP^*V^*，*=G、R、C。

2) 按如下严密公式进行方差分量估计：

$ {{\mathit{\boldsymbol{S}}}_{3 \times 3}}{{\mathit{\boldsymbol{\hat \theta }}}_{3 \times 1}} = {{\mathit{\boldsymbol{W}}}_{3 \times 1}} $

(11)

式中，

$ \begin{array}{l} \mathit{\boldsymbol{S}} = \\ \left[{\begin{array}{*{20}{c}} {{n_1}-2{\rm{tr}}({\mathit{\boldsymbol{N}}^{-1}}{\mathit{\boldsymbol{N}}^{\rm{G}}}) + {\rm{tr}}{{({\mathit{\boldsymbol{N}}^{-1}}{\mathit{\boldsymbol{N}}^{\rm{G}}})}^2}}&{{\rm{tr}}({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{G}}}{\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{R}}})}&{{\rm{tr}}({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{G}}}{\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{C}}})}\\ {{\rm{tr}}({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{G}}}{\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{R}}})}&{{n_2} - 2{\rm{tr}}({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{R}}}) + {\rm{tr}}{{({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{R}}})}^2}}&{{\rm{tr}}({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{R}}}{\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{C}}})}\\ {{\rm{tr}}({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{G}}}{\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{C}}})}&{{\rm{tr}}({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{R}}}{\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{C}}})}&{{n_3} - 2{\rm{tr}}({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{C}}}) + {\rm{tr}}{{({\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{N}}^{\rm{C}}})}^2}} \end{array}} \right] \end{array} $

$ {\mathit{\boldsymbol{\hat \theta }}} = \left[\begin{array}{l} \hat \sigma _{0, {\rm{G}}}^2\\ \hat \sigma _{0, {\rm{R}}}^2\\ \hat \sigma _{0, {\rm{C}}}^2 \end{array} \right], {\mathit{\boldsymbol{W}}} = \left[\begin{array}{l} {({{\mathit{\boldsymbol{V}}}^{\rm{G}}})^{\rm{T}}}{{\mathit{\boldsymbol{P}}}^{\rm{G}}}{{\mathit{\boldsymbol{V}}}^{\rm{G}}}\\ {({{\mathit{\boldsymbol{V}}}^{\rm{R}}})^{\rm{T}}}{{\mathit{\boldsymbol{P}}}^{\rm{R}}}{{\mathit{\boldsymbol{V}}}^{\rm{R}}}\\ {({{\mathit{\boldsymbol{V}}}^{\rm{C}}})^{\rm{T}}}{{\mathit{\boldsymbol{P}}}^{\rm{C}}}{{\mathit{\boldsymbol{V}}}^{\rm{C}}} \end{array} \right], {\mathit{\boldsymbol{N}}} = {{\mathit{\boldsymbol{N}}}^{\rm{G}}} + {{\mathit{\boldsymbol{N}}}^{\rm{R}}} + {\mathit{\boldsymbol{N}}^{\rm{C}}}, $

式中，n₁、n₂、n₃分别为GPS、GLONASS、BDS观测方程的个数。

3) 解得$\hat \sigma _{0, {\rm{G}}}^2$、$\hat \sigma _{0, \rm{R}}^2$、$\hat \sigma _{0, \rm{C}}^2$之后，重新调整权阵:

$ {({{\mathit{\boldsymbol{P}}}^*})_{k + 1}} = \frac{c}{{\hat \sigma _{0, *}^2}}{({{\mathit{\boldsymbol{P}}}^*})_k}, * = G, R, C $

(12)

式中，c为一常数，通常取$\hat \sigma _{0, *}^2$中的某个值。

4) 重复进行步骤1)~3)，直到$\hat \sigma _{0, {\rm{G}}}^2 \approx \hat \sigma _{0, {\rm{R}}}^2 \approx \hat \sigma _{0, {\rm{C}}}^2$为止，或通过必要的检验认为各类单位权方差之比约等于1为止。本文中，迭代的终止条件为:

$ \begin{array}{l} \max \left\{ {|\hat \sigma _{0, {\rm{G}}}^2-\hat \sigma _{0, {\rm{R}}}^2|, |\hat \sigma _{0, {\rm{G}}}^2-\hat \sigma _{0, {\rm{C}}}^2|, } \right.\\ \;\;\;\;\;\;\;\;\;\left. {|\hat \sigma _{0, {\rm{R}}}^2-\hat \sigma _{0, {\rm{C}}}^2|} \right\} \le 0.01 \end{array} $

(13)

4 算例及精度分析

采用GNSS三系统双频数据(GPS：L₁/L₂；GLONASS：L₁/L₂；BDS：B₁/B₂)，采样率为1 Hz，接收机为Propak6，天线类型为NOV703GGG，基准站为山东科技大学实验楼楼顶的S2ld，流动站为校园内点TaiY，其精确坐标均已知。数据采集时间为2015-05-26 02:00:00~03:59:59(UTC)。

4.1 多星座GNSS伪距观测值先验权比

多星座GNSS组合进行伪距单点定位时，采用Helmert方差分量估计方法计算GPS、GLONASS和BDS的伪距观测值方差的比值，结果见表 1。

为验证Helmert方差分量估计得到的伪距观测值权比的正确性，采用不同权比(GPS :GLONASS :BDS)进行组合单点定位。表 2为GPS/GLONASS/BDS组合伪距单点定位点位误差RMS统计。由表 2可以看出，组合伪距单点定位中，改变不同系统观测值的权比，定位精度有明显变化。GPS/GLONASS/BDS组合伪距单点定位时，其观测值的先验权比设置为5 :1 :1，组合定位精度最高，同时说明Helmert方差分量估计得到的GPS、GLONASS和BDS的伪距观测值权比是正确的。

4.2 多星座GNSS相位观测值先验权比

表 3为GPS/GLONASS/BDS组合基线解算中相位观测值的方差和权的比值。

为验证该权比的正确性，同样基于不同先验权比进行基线解算，点位误差结果见表 4。由表 4可以看出，GPS/GLONASS/BDS组合基线解算相位观测值的权比取不同值时，组合系统定位的固定解精度差异较小，该结论与GPS、GLONASS和BDS的相位观测值的精度相当一致。因此，GPS/GLONASS/BDS组合基线解算中，各系统相位观测值的权比可以简单设置为1 :1 :1。

5 结语

由于多GNSS系统间观测噪声、轨道精度的差异，采用经验的权比进行组合定位难以得到最优的结果。基于此，本文在GPS/GLONASS/BDS组合定位中引入Helmert方差分量估计，对GPS/GLONASS/BDS组合单点定位和基线解算中各系统观测值进行合理定权。为了验证该方法的正确性，基于同一组观测数据，采用不同的观测值权比进行组合定位验证，并统计N、E和U方向偏差的RMS和点位误差RMS来评价定位的性能。结果表明，使用Helmert方差分量估计确定的权比(本实验中伪距观测值的最优权比为5 :1 :1；相位观测值的最优权比为1 :1 :1)得到的多系统组合定位结果是最优的。在多频多模数据处理中，建议采用方差分量估计方法定权，可以有效提高多系统组合定位的精度和可靠性。