基于变分贝叶斯自适应鲁棒滤波的动对动相对定位算法

蔡 巍，陈明剑，邓 垦,3，王 洋，沈 洋，张擎天

（1.信息工程大学 地理空间信息学院，郑州 450001；2.西安测绘总站，西安 710054；3.77120 部队，成都 610000；4.河南北斗卫星导航平台有限公司，郑州 450001）

舰载机着舰、无人机编队飞行、装甲车编队行进等军事应用都涉及到动对动高精度相对定位技术，其特点是将基准站设置在移动的载体上，利用载波相位差分定位求解两个运动载体的相对位置。需要解决的关键问题是在动态情况下如何快速、准确地解算出整周模糊度[1]。在动对动定位条件下，观测条件变化多样，可见卫星数波动频繁，一旦发生周跳或卫星信号丢失，会导致模糊度固定率下降甚至是无法固定，需要重新搜索固定模糊度。因此，如何提高模糊度的正确性、可靠解和搜索效率是需要重点研究的问题。

北斗三号卫星导航系统（Beidou-3 Navigation Satellite System,BDS-3）具有播发五个频率观测值的功能，将多个频率观测值进行组合得到具有波长较长、电离层影响较弱以及低噪声等特性的组合观测值，可显著提高定位解算效率和模糊度固定率。伍劭实[2]利用北斗双频组合进行动对动相对定位，与单频解算相比，初始化时间用时更短。信息工程大学的许扬胤[3]围绕北斗三频信号在约束基线长度的情况下进行动对动相对定位技术研究。在短基线下给出了一种几何与无几何相结合的部分模糊度解算方案，模糊度的固定率可达到99.87%。刘濛濛[4]为了充分利用北斗三号B1I、B3I、B1C 和B2a 四频观测数据，采用无几何模型和几何模型相结合的方法依次求解两个超宽巷、三个宽巷和一个窄巷模糊度，并利用最小二乘模糊度降相关平差法（Least-square Ambiguity Decorrelation Adjustment,LAMBDA）方法搜索固定窄巷模糊度，最后使用模糊度固定的窄巷进行单历元定位。在短基线情况下，该方法可以有效进行超宽巷和宽巷模糊度的单历元固定，窄巷模糊度固定率高于99%，基本可以实现短基线单历元定位，并获得厘米级定位精度。王剑超[5]利用北斗三号B1I、B3I、B1C 和B2a 四频观测数据两两组合形成独立不相关的窄巷观测值联立，以位置参数和双差模糊度组成状态向量，结合卡尔曼滤波估计进行定位解算，在短基线情况下可以得到毫米级的定位结果，定位精度相较于北斗单频及三频方法在平面和高程方向均有不同程度的提高。随着观测频率的增加，可以组合出更多具有良好质量的观测值，较优的四频组合相对于传统的双频和三频的组合观测值具有波长较长、电离层放大因子较小、噪声放大因子较小、总噪声水平较小等特性；在短基线条件下，双差大气延迟误差对定位影响较小，且对流层延迟的大部分可以通过模型进行改正，因此可采用几何模型和无几何模型单历元固定整周模糊度。

以上所提出的动对动相对定位方法均是基于传统的双频和三频模糊度解算方法（Three Carrier Ambiguity Resolution,TCAR）求解模糊度，滤波方式均采用普通的卡尔曼滤波，该滤波方式在动态环境下不具有很好的自适应性和鲁棒性。四频组合观测值相较于传统的双频和三频组合观测值在模糊度固定率、定位精度方面都有提高，而将四频组合观测值用于动对动相对定位的方法较少。四频数据组合可构造更多的弱电离层延迟和小观测噪声影响的组合观测值，对于提高模糊度固定率有重要意义[6-8]，因此，本文采用“宽巷+窄巷+伪距”的四频信号无几何模糊度解算模型，该模型提高了传统的无几何模型模糊度解算效率和可靠性。为适应动态环境的干扰，本文提出变分贝叶斯自适应鲁棒滤波算法，兼顾滤波过程中的自适应性与鲁棒性。最后通过实验分析动对动定位方法解算模糊度固定率和基线解算精度。

1 BDS 的四频线性组合

利用BDS-3 的四个频率观测数据，通过不同频率载波相位观测值间的线性组合，可以获得不同波长、不同电离层延迟误差放大系数、不同观测噪声放大因子等组合观测值。本文采用的BDS-3 的四个频率分别是B1I、B1C、B3I 和B2a，组合后的频率、波长和模糊度表达式为：

其中，ϕ1、ϕ2、ϕ3、ϕ4分别为四个频率信号的双差相位；P1、P2、P3、P4分别为四个频率信号的伪距观测值；a1、b2、c3、d4分别为伪距的组合系数；分别为组合后的双差伪距和相位观测值，其双差观测方程为：

将伪距组合系数代入式(6)即可得到伪距组合观测值的电离层延迟误差放大系数。

假设各个频点的观测值噪声都是相等且相互独立的，则对应的多频组合观测噪声表达式为：

将伪距组合系数代入式(8)即可得到伪距组合观测噪声放大因子。

在短基线条件下，双差大气延迟误差对定位影响较小，且对流层延迟的大部分可以通过模型进行改正，将式(5)减去式(4)就可得到基于四频组合观测值的无几何观测方程，其表达式为：

可利用无几何观测方程求解超宽巷和宽巷的模糊度。

2 BDS 常用的四频线性组合

随着不同频点观测值的增加，可以利用各个频点观测值进行线性组合，得到更多质量好的组合观测值。在实际应用中，通常选择波长较长、大气延迟和观测噪声较小的组合观测值。按照文献[9-10]对不同波长的观测值进行分类，可分为三类：1）超宽巷λ≥2.93 m ；2）宽巷0.75 m ≤λ＜2.93m ；3）窄巷0.1m ≤λ＜0.75m。不同组合系数在不同长度基线下所表现出的噪声水平是不同的，在选取最佳组合系数时，引入总噪声水平对其进行评价，其表达式为：

其中，δ为观测总噪声水平；δorb为轨道误差；δtrop为对流层延迟误差；δI1与δI2分别为一阶和二阶电离层延迟误差；δε为载波相位观测噪声。当基线长度为中短基线时（长度≤100 km），以上各个误差的取值分别为：0.005 m、0.01 m、0.1 m、0.005 m、0.005 m。表1 按照波长、电离层延迟误差放大系数、观测噪声放大因子和总噪声水平四个方面选出最佳的超宽巷、宽巷和窄巷的组合系数。

表1 BDS-3 四频系数组合Tab.1 Coefficient combinations of BDS-3 quad-frequency observations

通过表1 可知，（0,0,1,-1）为超宽巷最佳组合系数，（0,1,-3,2）、（1,0,-2,1）和（0,1,-1,0）为最佳的宽巷组合系数，（2,2,0,-3）为最佳的窄巷组合系数。

3 BDS 四频组合观测值动对动相对定位算法

3.1 基于GF 模型单历元超宽巷/宽巷模糊度固定

通过式(9)可知双差伪距观测值与双差载波相位观测值之差可得超宽巷模糊度的浮点解，只有选取最优的伪距组合观测值，其超宽巷模糊度浮点解的精度才最优。在伪距观测噪声为0.5 m 的条件下，以无几何（Geometric-Free,GF）模型的总噪声水平为标准选取最优伪距组合观测值，通过比较，GF 模型Δ∇P(1,1,1,1)-Δ∇ϕ(0,0,1,-1)的总噪声水平最低，总噪声水平为0.0834 周，远远小于0.5 周。理论上利用GF 组合模型求出的超宽巷模糊度通过四舍五入直接取整就可以进行单历元固定。

通过Δ∇P(1,1,1,1)-Δ∇ϕ(0,0,1,-1)组合可求出超宽巷模糊度浮点解Δ∇N(0,0,1,-1)，对其进行四舍五入取整可得超宽巷模糊度整数解。本文以采集的47 km 长基线数据为例，两颗卫星C35-C22 用GF 模型双差后得到的浮点解Δ∇N(0,0,1,-1).进行相邻历元间差分，其结果见图1。

图1 超宽巷（0,0,1,-1）相邻历元间模糊度差值Fig.1 Ultra-wide lane (0,0,1,-1) ambiguity difference between adjacent epoch

通过图1 可以看出，相邻历元间超宽巷模糊度浮点解的差值远小于±0.5 周，可直接采用四舍五入单历元取整的方法，可靠地固定超宽巷模糊度。

图2 宽巷（0,1,-3,2）相邻历元间模糊度差值Fig.2 Wide lane (0,1,-3,2) ambiguity difference between adjacent epoch

图3 宽巷（1,0,-2,1）相邻历元间模糊度差值Fig.3 Wide lane (1,0,-2,1) ambiguity difference between adjacent epoch

图4 宽巷（0,1,-1,0）相邻历元间模糊度差值Fig.4 Wide lane (0,1,-1,0) ambiguity difference between adjacent epoch

通过图2-4 可以看出，利用模糊度固定的超宽巷观测值作为约束条件，得到相邻历元间宽巷模糊度浮点解的差值绝对值小于0.5 周，可直接采用四舍五入单历元取整方法，可靠地固定宽巷模糊度。

3.2 变分贝叶斯自适应鲁棒滤波模型

卡尔曼滤波在噪声参数已知并保持不变的情况下解算性能优越[11]，然而在实际环境中，将噪声设置为固定的高斯白噪声具有一定的局限性。当目标处于运动状态下，系统受到外部环境未知噪声的干扰，导致观测噪声在滤波过程中是变化的，此时卡尔曼滤波解算性能受到影响，并不能体现其优越性。将变分贝叶斯估计与卡尔曼滤波相结合，对量测噪声方差进行多次迭代估计，让量测噪声方差适应动态环境下的滤波。

滤波器的状态及量测方程为：

其中，Xk为当前时刻状态向量；k为当前历元时刻，k-1为上一个历元时刻；Φk,k-1为两个时刻状态转移矩阵；ωk-1为系统噪声矩阵；Hk为线性化观测矩阵；vk表示测量噪声；Zk为量测向量，其表达式为：

其中，ι=tk-tk-1表示接收机的采用间隔；有m-1个载波相位双差方程。

Hk表达式为：

其中，λ为组合观测值波长；E表示三维坐标改正值的系数阵，其表达式为：

其中，e1表示参考星单位视线向量；e2-ej(j=2,…,m)为非参考星单位视线向量。

Qk是ωk的协方差矩阵，在动态下的卡尔曼滤波表达式为：

状态向量Xk包括：

对于先验观测误差方差阵中第i个方差Rk,i，计算公式为：

其中，eli代表对应第i颗卫星的高度角；当观测值为窄巷组合观测值时，a和b为载波相位观测误差；当观测值为宽巷组合观测值时，a和b为伪距观测误差；ε为基线向量观测误差；d为卫星钟稳定误差。

变分贝叶斯滤波的目的是计算出后验分布p(Xk,Rk|Z1:k)，选择系统状态量Xk和观测误差方差Rk作为待估对象。迭代估计的步骤可以分为如下几步：

（1）初始化：先从先验分布开始，先验分布为p(Xk-1,Rk-1|Z1:k-1)，其中Xk-1、Z1:k-1均为初始值，Rk-1为先验观测误差方差阵。

（2）预测：对状态和方差进行预测。

（3）量测更新：根据贝叶斯公式，可以得到后验分布：

由于Xk服从高斯分布，其共轭分布为逆伽马分布，文献[12]也是将观测误差的方差分布选择为逆伽马分布。但这种方法在更新阶段中，状态和观测误差的方差是耦合的，需要用一个新的分布来近似真实的后验分布：

其中，QX(Xk)和QR(Rk)分别是基于Xk和Rk的概率密度分布，其表达式为：

变分贝叶斯方法通过最小化两个分布间的Kullback-Leibler（KL）散度，使构造的变分分布逼近真实的联合后验分布。通过最小化真实后验分布与近似后验分布的KL 散度距离，QX(Xk)和QR(Rk)近似表达式为[13-15]：

其中，mk为待估参数，是近似得到的新分布的状态量；Pk为状态协方差矩阵为观测误差方差阵Rk的对角线元素；Inv-Gamma(·)表示逆伽玛分布，αk,i和βk,i为逆伽玛分布的待估参数。由式(23)(24)可看出，状态Xk服从高斯分布，观测误差方差Rk服从逆伽马分布。

通过对比基于QX(Xk)和QR(Rk)的期望值得出参数估计为：

对观测误差方差的鲁棒化采用M 估计的方式。

首先由量测方程Zk-HkXk=vk，得到量测噪声vk。然后标准化量测噪声：

其中，表示标准化后的量测噪声。

构造权重函数Ψk,i，其表达式为：

其中，γ为调节因子，是代价函数的门限，体现了对外部干扰情况的处理程度，一般取1.345。M 估计的鲁棒化思想是重加权平均，其方法是将权矩阵Ψ进行分块处理，该权阵的元素由Ψk,i构成，分为状态量权阵Ψx和观测权阵Ψy，其表达式为[16]：

则观测误差方差阵经过鲁棒化得到的表达式为：

在参数估计中一般循环迭代三次，每个历元通过构造的滤波器得到坐标参数和模糊度的浮点解，然后利用LAMBDA 方法[17]对窄巷模糊度进行单历元固定，再将窄巷模糊度固定解回代窄巷组合观测方程，最后得到位置坐标的固定解。

4 实验分析

为验证本文动对动相对定位算法的可行性，采用一条河南省内 CORS 站网的基线，基线长度为47.2 km，观测日期为2022 年9 月1 日8:00-9:00，采样间隔为1 s。解算过程中，设定卫星截止高度角为20 °，模糊度的Ratio 阀值采用3.0。

实验一：采用静态数据模拟动态数据进行，在静态数据中加入随机游走噪声。基准站位置采用单点定位，模糊度固定方式采用fix and hold 方式。实验所用到的状态先验方差和过程噪声见表2。由于本实验只是验证算法的可行性，因此不加入其他算法进行比较，图5 给出两个测站之间满足解算条件的共视卫星数量，图6-8 给出窄巷组合定位结果与真值之间的偏差，其结果统计见表3。

图5 实验一满足解算条件的共视卫星数Fig.5 The number of common-view satellites satisfying the solution condition in the experiment 1

图6 实验一解算结果在E 方向与真值的偏差Fig.6 The deviation of the result from the true value in the E direction in the experiment 1

图7 实验一解算结果在N 方向与真值的偏差Fig.7 The deviation of the result from the true value in the N direction in the experiment 1

图8 实验一解算结果在U 方向与真值的偏差Fig.8 The deviation of the result from the true value in the U direction in the experiment 1

表2 状态先验方差及过程噪声Tab.2 State prior variance and process noise

表3 实验一解算结果Tab.3 Solution results in the experiment 1

从图6-8 和表3 可以看出，采用本文算法在静态数据模拟动态数据的情况下，基线的E、N 和U 方向定位结果与真值的偏差均能达到厘米级，模糊度固定率可以达到99%以上。

实验二：为验证本文算法在真实的动态环境下的解算性能，动态实验于2023 年6 月5 日在信息工程大学的两个操场间进行，实验时间为 30 min左右。实验设备为3 台司南接收机，将一台司南接收机设置为固定基准站，并实时通过网络RTK 获取定位数据，与其他两台用于移动的接收机构成静动型基线，将两台移动接收机精确的定位结果作为真值与动对动实验结果比较。两台移动的接收机分别绑在两辆电动车上，在两个不同的操场上相对运动，基线长度为0～400 m，将其中一台接收机充当移动的基准站，其定位结果采用单点定位，另外一台接收机充当移动的流动站，与其构成动态条件下的基线，定位结果采用动对动相对定位算法。采用四频普通卡尔曼滤波定位解算结果和RTKLIB 三频非组合动基线解算结果作为其比较对象，实验参数设置与实验一相同，运动轨迹见图9。

图9 动对动实验轨迹Fig.9 Motion to motion experiment trajectory

图10 给出了实验二的两个测站之间满足解算条件的共视卫星数量，图11-14 给出分别采用三种不同算法解算定位结果分量、基线长度与真值之间的偏差，其对比结果见表4。

图10 实验二满足解算条件的共视卫星数Fig.10 The number of common-view satellites satisfying the solution condition in the experiment 2

图11 实验二解算结果在E 方向与真值的偏差Fig.11 The deviation of the result from the true value in the E direction in the experiment 2

图12 实验二解算结果在N 方向与真值的偏差Fig.12 The deviation of the result from the true value in the N direction in the experiment 2

图13 实验二解算结果在U 方向与真值的偏差Fig.13 The deviation of the result from the true value in the U direction in the experiment 2

图14 实验二解算基线长度与真值的偏差Fig.14 The deviation of the baseline length from the true value in the experiment 2

表4 实验二不同算法解算结果比较Tab.4 Comparison of the results of different algorithms in the experiment 2

从图11-14 和表4 可以看出，在动态环境下，采用北斗四频的变分贝叶斯自适应鲁棒滤波算法，在东、北、天三个方向的定位精度均优于2 cm，基线长度与真值偏差达到毫米级，模糊度固定率达到100%，证明采用四频组合观测值的贝叶斯滤波要优于四频组合观测值的卡尔曼滤波；随着不同频点观测值的增加，可以利用各个频点观测值进行线性组合，得到更多质量好的组合观测值，有利于模糊度的固定。

5 结论

本文采用BDS-3 的B1I、B1C、B3I 和B2a 四频信号进行组合，得到不同的组合观测值，再结合变分贝叶斯自适应鲁棒滤波算法，进行动对动相对定位。对可见卫星数、模糊度固定率和定位精度进行了统计。四频组合相对于三频可以得到波长更长、观测噪声更小的组合观测值，进一步提高了组合观测值模糊度的固定效率。在静态模拟动态实验的3600 个历元中，模糊度的固定率达到了99.7%，定位精度达到厘米级；在更加复杂的动态环境下，模糊度的固定率达到了100%，定位结果与真值的偏差优于2 cm，验证了基于北斗四频的变分贝叶斯自适应鲁棒滤波算法可以满足高精度的动对动相对定位。