BDS/GPS紧组合相比松组合能够提高导航定位的可用性、可靠性和精度，但在使用紧组合模型时必须考虑接收机端差分系统间偏差(DISB)的影响^[1]。Paziewski等^[2]、Odijk等^[3-4]研究了GPS L₁/Galileo E₁紧组合模型；张小红等^[1]、Wu等^[5]、周英东等^[6]研究了BDS B₂/Galileo E_5b紧组合模型，利用最小二乘估计紧组合中的DISB参数，并分析DISB的特性。此外，Odolinski等^[7]研究了BDS、Galileo、QZSS和GPS 4系统紧组合单频RTK，但由于BDS B₁与其他系统的频率不一致，文中并未详细给出BDS与其他系统的紧组合算法以及DISB的估计结果。

相比GPS L₁/Galileo E₁和BDS B₂/Galileo E_5b的紧组合来说，BDS B₁/GPS L₁紧组合由于存在频率不一致问题，估计的载波DISB小数部分易受到参考卫星站间单差模糊度的影响。楼益栋等^[8]在研究BDS B₁/GPS L₁紧组合时，利用伪距计算的站间单差模糊度经多历元平滑后，改正由频率不一致引入的参考卫星站间单差模糊度，然后再由最小二乘估计DISB参数；隋心等^[9]提出利用卡尔曼滤波算法，同时估计站间单差模糊度和DISB参数。但这2种算法都只估计了载波DISB小数部分，没有考虑载波DISB整数部分。然而由于BDS和GPS的信号频率不一致，在估计载波DISB小数部分时，其整数部分会影响参考卫星站间单差模糊度的精度。另外，Tian等^[10]针对GPS L₁/Galileo E₁紧组合中载波DISB的估计问题，提出一种粒子滤波估计算法。本文针对BDS B₁/GPS L₁紧组合推导了其相对定位模型，考虑到BDS B₁/GPS L₁紧组合中频率不一致的问题，提出联合卡尔曼滤波与粒子滤波算法估计BDS B₁/GPS L₁紧组合中的DISB，并通过实验对比分析卡尔曼滤波算法和本文算法的效果。

1 BDS/GPS紧组合相对定位模型

在相对定位时，假设接收机r(=1, 2)能同时接收到BDS和GPS卫星的伪距和载波，则对于任意卫星p和接收机r，可得非差观测方程：

$ \begin{array}{*{20}{c}} {P_r^p = \rho _r^p + c\left( {{\rm{d}}{T_r} - {\rm{d}}{t^p}} \right) + }\\ {{\rm{d}}{r_r} - {\rm{d}}{s^p} + I_r^p + T_r^p + {\varepsilon _p}} \end{array} $

(1)

$ \begin{array}{*{20}{c}} {{\lambda ^p}\varphi _r^p = \rho _r^p + c\left( {{\rm{d}}{T_r} - {\rm{d}}{t^p}} \right) + }\\ {{\rm{b}}{{\rm{r}}_r} - {\rm{b}}{{\rm{s}}^p} + {\lambda ^p}N_r^p - I_r^p + T_r^p + {\varepsilon _\varphi }} \end{array} $

(2)

式中，P为伪距观测值，ρ为卫星到接收机的几何距离，c为光在真空中的传播速度，dT为接收机钟差，dt为卫星钟差，dr为接收机对伪距的硬件延迟，ds为卫星对伪距的硬件延迟，I为电离层延迟，T为对流层延迟，ε_p为伪距测量噪声，λ为载波波长，φ为载波观测值，br为接收机对载波的硬件延迟，bs为卫星对载波的硬件延迟，N为整周模糊度，ε_φ为载波观测噪声。

在接收机间求差可得站间单差观测方程：

$ \Delta P_{12}^p = \Delta \rho _{12}^p + c\Delta {\rm{d}}{T_{12}} + \Delta {\rm{d}}r_{12}^p + \Delta \varepsilon _p^p $

(3)

$ \begin{array}{*{20}{c}} {{\lambda ^p}\Delta \varphi _{12}^p = \Delta \rho _{12}^p + c\Delta {\rm{d}}{T_{12}} + }\\ {\Delta {\rm{br}}_{12}^p + {\lambda ^p}\Delta N_{12}^p + \Delta \varepsilon _\varphi ^p} \end{array} $

(4)

式中，Δ(·)₁₂^p为站间单差运算符号，上标p=G、B分别为GPS和BDS卫星。站间单差可以消除与卫星相关的卫星钟差和卫星端硬件延迟等，削弱电离层延迟和对流层延迟误差等大气延迟；对于短基线，残留的大气延迟可以忽略不计。另外还需指出，当基线两端接收机类型不同时，BDS的GEO卫星与IGSO/MEO卫星之间还可能存在ISTB(inter-satellite type bias)的影响，在数据处理时需提前改正^[11-12]。

与常规的BDS/GPS松组合相对定位模型不同，紧组合模型在2个系统中只选取了1颗参考卫星。假设选取GPS中的G₁卫星作为参考卫星，则可组成站星双差观测方程：

$ \Delta \nabla P_{12}^{{G_1}{G_i}} = \Delta \nabla \rho _{12}^{{G_1}{G_i}} + \Delta \nabla \varepsilon _p^{GG} $

(5)

$ \begin{array}{*{20}{c}} {{\lambda ^G}\Delta \nabla \varphi _{12}^{{G_1}{G_i}} = }\\ {\Delta \nabla \rho _{12}^{{G_1}{G_i}} + {\lambda ^G}\Delta \nabla N_{12}^{{G_1}{G_i}} + \Delta \nabla \varepsilon _p^{GG}} \end{array} $

(6)

$ \Delta \nabla P_{12}^{{G_1}{B_j}} = \Delta \nabla \rho _{12}^{{G_1}{B_j}} + \Delta {\rm{d}}r_{12}^{GB} + \Delta \nabla \varepsilon _p^{GB} $

(7)

$ \begin{array}{*{20}{c}} {{\lambda ^B}\Delta \varphi _{12}^{{B_j}} - {\lambda ^G}\Delta \varphi _{12}^{{G_1}} = \Delta \nabla \rho _{12}^{{G_1}{B_j}} + }\\ {\Delta \nabla {\rm{br}}_{12}^{GB} + \left( {{\lambda ^B}\Delta N_{12}^{{B_j}} - {\lambda ^G}\Delta N_{12}^{{G_1}}} \right) + \Delta \nabla \varepsilon _\varphi ^{GB}} \end{array} $

(8)

式中，Δ∇(·)₁₂^pq为卫星p、q和接收机1、2间的双差运算符号，G_i=G₁, G₂…G_n为GPS观测卫星，B_j=B₁, B₂…B_m为BDS观测卫星，Δdr₁₂^GB为伪距DISB，Δ∇br₁₂^GB为载波DISB。DISB是在系统间双差观测方程中，由于2颗卫星所属的系统不同而在接收机端引起的系统间偏差。

在估计出DISB小数部分后，可用估计的DISB改正紧组合，当采用卡尔曼滤波估计站间单差模糊度的方式处理改正DISB后的紧组合相对定位模型时，方程可表示为：

$ \begin{array}{*{20}{c}} {{\lambda ^B}\Delta \varphi _{12}^{{B_j}} - {\lambda ^G}\Delta \varphi _{12}^{{G_1}} - {\lambda ^B}{F_{{\rm{br}}}} = }\\ {\Delta \nabla \rho _{12}^{{G_1}{B_j}} + {\lambda ^B}{N_{{\rm{br}}}} + {\lambda ^B}\nabla N_{12}^{{B_j}} - {\lambda ^G}\Delta N_{12}^{{G_1}} + }\\ {\Delta \nabla \varepsilon _p^{GB} = \Delta \nabla \rho _{12}^{{G_1}{B_j}} + {\lambda ^B}\left( {\nabla N_{12}^{{B_j}} + N_{{\rm{br}}}^B} \right) - }\\ {{\lambda ^G}\left( {\nabla N_{12}^{{G_1}} - N_{{\rm{br}}}^G} \right) + \Delta \nabla \varepsilon _p^{GB} = \Delta \nabla \rho _{12}^{{G_1}{B_j}} + {\lambda ^B} \cdot }\\ {\left( {\Delta \nabla N_{12}^{{G_1}{B_j}} + {N_{{\rm{br}}}}} \right) + \left( {{\lambda ^B} - {\lambda ^G}} \right) \cdot }\\ {\left( {\nabla N_{12}^{{G_1}} - N_{{\rm{br}}}^G} \right) + \Delta \nabla \varepsilon _p^{GB}} \end{array} $

(9)

式中，F_br为已知的载波DISB小数部分，N_br为载波DISB整数部分，N_br^B和N_br^G为载波DISB整数部分对站间单差模糊度的影响，且λ^BN_br=λ^BN_br^B+λ^GN_br^G。式(9)中，最后一个等号表示由模糊度浮点解到固定双差模糊度，由于双差模糊度被固定为整数，则等号右边会多出载波DISB整数部分对参考卫星站间单差模糊度的影响，故方程中还会残留其整数部分与频率不一致引起的小数偏差，因此在估计载波DISB小数部分时还需提前考虑该偏差。

2 DISB估计算法 2.1 卡尔曼滤波估计算法

载波DISB和频率不一致参数的存在，导致方程秩亏，并且载波DISB与模糊度参数紧密联系在一起，很难分离。因此，利用卡尔曼滤波算法估计DISB参数时，首先要对方程进行参数重整：

$ \begin{array}{*{20}{c}} {{\lambda ^B}\Delta \varphi _{12}^{{B_j}} - {\lambda ^G}\Delta \varphi _{12}^{{G_1}} = \Delta \nabla \rho _{12}^{{G_1}{B_j}} + }\\ {\left( {\Delta \nabla {\rm{br}}_{12}^{GB} + {\lambda ^B}\Delta \nabla N_{12}^{{G_1}{B_1}}} \right) + {\lambda ^B} \cdot }\\ {\Delta \nabla N_{12}^{{B_1}{B_j}} + \left( {{\lambda ^B} - {\lambda ^G}} \right)\Delta N_{12}^{{G_1}} + \Delta \nabla \varepsilon _p^{GB}} \end{array} $

(10)

考虑到存在频率不一致引起的ΔN₁₂^G₁参数，故将Δ∇N₁₂^G₁G_i和Δ∇N₁₂^B₁B_j均以站间单差模糊度的方式进行估计，且在方程中将Δ∇br₁₂^GB+λ^B·Δ∇N₁₂^G₁B₁当作一个λ^BΔ∇$\widetilde N$₁₂^G₁B₁参数，然后一起由卡尔曼滤波估计。在卡尔曼滤波估计出站间单差模糊度浮点解后，将站间单差模糊度变换为双差模糊度^[13]，并用LAMBDA算法^[14]对双差模糊度Δ∇N₁₂^G₁G_i和Δ∇N₁₂^B₁B_j进行固定，计算固定模糊度后的坐标和DISB参数；在估计出Δ∇$\widetilde N$₁₂^G₁B₁后，对其直接取整后的小数部分即为载波DISB的小数部分，而载波DISB整数部分则可被双差模糊度吸收。另外，伪距DISB可以和坐标参数、模糊度参数一起估计。该算法同时解决了频率不一致以及DISB与模糊度强相关的问题，但并没有提前改正载波DISB整数部分与频率差异产生的小数偏差。

2.2 联合卡尔曼滤波与粒子滤波估计算法

本文联合算法的基本思想是：首先基于卡尔曼滤波和粒子滤波算法估计出载波DISB的小数部分，然后再采用卡尔曼滤波算法一同估计坐标参数、伪距DISB参数以及模糊度参数。由于粒子滤波只需估计载波DISB不足1周的部分，故载波DISB的取值范围可假设为[-0.5, 0.5]周，约为[-0.1, 0.1] m。给定采样间隔为1 mm，可以得到N=200个粒子x₀ⁱ，给定每个粒子的初始权重为w₀ⁱ=1/N，其中i=1, 2, …, N。

当载波DISB小数部分已知时，粒子滤波的观测方程为：

$ \mathit{\boldsymbol{Z}} = \mathit{\boldsymbol{HX}} + \mathit{\boldsymbol{v}} $

(11)

式中，Z为观测值减去计算值残差，H为设计矩阵，X为状态向量，v为测量噪声，且有X=[x y z Δ∇dr ΔN₁₂^G₁ … ΔN₁₂^G_n ΔN₁₂^B₁ … ΔN₁₂^B_m]^T。

参考文献[10, 15]中的粒子滤波算法，先根据概率函数更新每个粒子的权重。对于概率函数的选取，文献[10, 15]通过实验证明，当一个越接近正确值的DISB值代入紧组合模型中时，通过LAMBDA算法固定双差模糊度时的Ratio值越大，因此可根据Ratio值的大小来判断DISB的正确性，故本文采用双差模糊度固定时的Ratio值来构造观测值的概率函数。将每个粒子代入方程中，由卡尔曼滤波计算出站间单差模糊度浮点解，并将站间单差模糊度投影变换为双差模糊度，再用LAMBDA算法固定双差模糊度，计算出对应的Ratio值。于是可得每个粒子的概率函数：

$ p\left( {{b_k}\left| {\mathit{\boldsymbol{x}}_k^i} \right.} \right) = \frac{{{\rm{Rati}}{{\rm{o}}^i}}}{{\sum\limits_{j = 1}^N {{\rm{Rati}}{{\rm{o}}^j}} }} $

(12)

需要注意的是，每个粒子代入卡尔曼滤波计算出对应的Ratio值即可，此时状态参数不需要更新。根据概率函数更新每个粒子的权重并归一化：

$ \hat w_k^i = \frac{{w_{k - 1}^ip\left( {{b_k}\left| {\mathit{\boldsymbol{x}}_k^i} \right.} \right)}}{{\sum\limits_{j = 1}^N {w_{k - 1}^jp\left( {{b_k}\left| {\mathit{\boldsymbol{x}}_k^j} \right.} \right)} }} $

(13)

更新状态参数及其方差：

$ {{\mathit{\boldsymbol{\hat x}}}_k} \approx \sum\limits_{i = 1}^N {\hat w_k^i\mathit{\boldsymbol{x}}_k^i} $

(14)

$ {{\hat p}_k} \approx \sum\limits_{i = 1}^N {\left( {\mathit{\boldsymbol{x}}_k^i - {{\mathit{\boldsymbol{\hat x}}}_k}} \right){{\left( {\mathit{\boldsymbol{x}}_k^i - {{\mathit{\boldsymbol{\hat x}}}_k}} \right)}^{\rm{T}}}\hat w_k^i} $

(15)

每个历元计算出载波DISB后，将载波DISB估值${\mathit{\boldsymbol{\hat x}}_k}$代入观测方程中，然后再用卡尔曼滤波一同估计坐标参数、站间单差模糊度参数和伪距DISB参数，并将状态参数更新到下一个历元。其中，卡尔曼滤波状态方程可表示为：

$ {\mathit{\boldsymbol{X}}_k} = {\mathit{\boldsymbol{F}}_{k\left| {k - 1} \right.}}{\mathit{\boldsymbol{X}}_{k - 1}} + \mathit{\boldsymbol{w}} $

(16)

式中，X_k和X_k-1为k和k-1时刻的状态向量，F_k|k-1为状态转移矩阵，w为预测噪声。

一次粒子滤波计算完成后，每个粒子的验后权不再是相等的，靠近正确值附近的粒子权重较大，远离正确值的粒子权重较小。经多次计算后，权重较大的有效粒子数越来越少，因此，需要对粒子进行重采样^{[10, 15]}。计算过程如下：

构造一个随机序列：

$ \left\{ {{u_i} = \frac{{\left( {i - 1} \right) + {{\tilde u}_i}}}{N}} \right\}_{i = 1}^N $

(17)

式中，${\tilde u_i}$为服从U(0, 1)的随机分布。

计算累积权序列：

$ \left\{ {\tilde w_k^j = \sum\limits_{h = 1}^j {\hat w_k^h} } \right\}_{j = 1}^N $

(18)

比较$\tilde \omega _k^j$与u_i的大小，提取当$\tilde \omega _k^j > {u_i}$时的粒子，并将新获得的N个粒子重新给定权w_kⁱ=1/N。

需要说明的是，文献[10, 15]提出的重采样步骤计算量较大，当所有粒子的权重满足$1/\sum\nolimits_{i = 1}^N {\left( {{{\hat \omega }^i}} \right)} \ge 2/3N$时，可以不进行重采样。

由于紧组合中的DISB参数在短期内是稳定不变的，因此粒子滤波的预测方程可表示为：

$ {\mathit{\boldsymbol{x}}_k} = {\mathit{\boldsymbol{x}}_{k - 1}} + \mathit{\boldsymbol{w}} $

(19)

式中，x为载波DISB粒子，w为预测噪声。下一个历元的粒子即可通过方程进行预测。

3 实验分析

选取2016年年积日001天澳大利亚科廷大学的4条短基线数据进行实验，观测时长为24 h，采样间隔为30 s，详细信息见表 1。

实验中选取高度角最高的GPS卫星作为参考卫星，并设置卫星截止高度角为10°，然后利用2种算法估计BDS/GPS紧组合中的DISB。图 1给出了基线两端接收机类型和固件版本都相同的基线1和2由2种算法估计的DISB结果；图 2给出了基线两端接收机类型不同的基线3和4的DISB结果。由上节的估计算法可知，2种算法估计伪距DISB都是采用卡尔曼滤波，故其估计的伪距DISB基本一致，因此图 1和图 2中仅给出了卡尔曼滤波算法估计的伪距DISB结果。

2种算法的估计结果表明，BDS和GPS之间的DISB是稳定的，因此在BDS/GPS紧组合时可以提前校正DISB参数。此外，在卡尔曼滤波估计算法中，载波DISB是由估计的Δ∇$\tilde N$₁₂^G₁B₁分离而得，当G₁、B₁这2颗卫星有周跳发生时，需重新初始化Δ∇$\tilde N$₁₂^G₁B₁参数，这会导致估计的载波DISB不稳定，如图 1基线1和图 2基线3在大约3 h处明显波动；而本文算法却不会受到周跳的影响，稳定性更好。

另一方面，比较图 1和图 2中2种算法的DISB可以发现，无论基线两端接收机类型是否相同，2种算法估计的载波DISB小数部分都可能会存在一定偏差。从上节的理论分析可以看到，这一偏差主要是由载波DISB整数部分与频率差异引起的小数偏差。在卡尔曼滤波估计载波DISB小数部分时，其整数部分可被Δ∇N₁₂^G₁B_j直接吸收，此时，

$ \begin{array}{*{20}{c}} {{\lambda ^B}\Delta \varphi _{12}^{{B_j}} - {\lambda ^G}\Delta \varphi _{12}^{{G_1}} = \Delta \nabla \rho _{12}^{{G_1}{B_j}} + {\lambda ^B} \cdot }\\ {{{\tilde F}_{{\rm{br}}}} + {\lambda ^B}\left( {\Delta \nabla N_{12}^{{G_1}{B_1}} + {N_{{\rm{br}}}}} \right) + {\lambda ^B} \cdot }\\ {\Delta \nabla N_{12}^{{B_1}{B_j}} + \left( {{\lambda ^B} - {\lambda ^G}} \right)\Delta N_{12}^{{G_1}} + \Delta \nabla \varepsilon _\varphi ^{GB}} \end{array} $

(20)

而联合算法在估计DISB小数部分时，其整数部分最后虽然可被Δ∇N₁₂^G₁B_j吸收，但是在估计过程中，其整数部分会影响到参考卫星站间单差模糊度。比较式(20)和式(9)可得：

$ {{\tilde F}_{{\rm{br}}}} = {F_{{\rm{br}}}} - \left( {{\lambda ^B} - {\lambda ^G}} \right)N_{{\rm{br}}}^G $

(21)

式中，$\tilde F$_br为卡尔曼滤波估计的载波DISB小数部分，F_br为本文算法估计的载波DISB小数部分。由式(21)可知，当载波DISB存在整数部分时，2种算法估计的结果会存在由载波DISB整数部分和频率不一致共同引起的偏差。为验证这一理论推导结果，实验给出了4条基线在采用2种算法估计DISB时的参考卫星站间单差模糊度差值序列，即式(21)中N_br^G的时间序列，如图 3所示。结合图 1和图 2的2种算法估计的载波DISB小数部分均值，并考虑估计精度可以发现，4条基线结果均满足等式(21)。这表明载波DISB的整数部分会影响参考卫星站间单差模糊度的估计，从而影响载波DISB小数部分的估计结果。

为进一步验证2种算法估计结果的准确性和可靠性，将2种算法估计的DISB均值代入BDS/GPS紧组合相对定位模型中，分别设置卫星截止高度角为10°和30°，然后采用动态单历元模式解算这4条基线数据，并设定Ratio值阈值为3，统计固定解的RMS、平均Ratio值以及模糊度固定率，详细情况如表 2和表 3所示。

由表 2可知，在卫星截止高度角为10°时，由于观测卫星数较多，2种算法估计的DISB差异对紧组合定位结果影响不明显。由表 3可知，在卫星截止高度角为30°时，2种算法估计的载波DISB小数部分的差异对定位结果有比较明显的影响，且用本文算法估计的载波DISB小数部分改正紧组合后的定位精度，其模糊度固定率和可靠性均优于卡尔曼滤波。其中，对于差异最大的基线4来说，用本文算法估计的结果改正紧组合后的固定解RMS提高了10%~20%。

4 结语

本文推导了BDS B₁/GPS L₁紧组合相对定位模型，详细研究了BDS B₁/GPS L₁紧组合中DISB参数估计问题。考虑到BDS B₁/GPS L₁紧组合中存在频率不一致的问题，在估计载波DISB小数部分时还需考虑其整数部分与频率差异产生的小数偏差，本文提出联合卡尔曼滤波与粒子滤波算法估计BDS B₁/GPS L₁紧组合中的DISB，并利用4条实测短基线数据对比分析本文算法和卡尔曼滤波算法估计DISB的性能。结果表明，载波DISB的整数部分会影响参考卫星站间单差模糊度的估计，从而影响DISB小数部分的估计结果。将2种算法估计的DISB代入紧组合模型中进行单历元相对定位表明，相比用卡尔曼滤波估计结果改正DISB，用本文算法估计结果改正DISB后的紧组合在模糊度固定率、模糊度固定可靠性以及固定解RMS方面均有一定提高，尤其是当载波DISB整数部分较大时，固定解精度可提高10%~20%。

致谢: 感谢科廷大学GNSS中心(Curtin GNSS Research Centre)提供的短基线数据。