﻿ 城市环境下BDS/GPS单频RTK定位算法研究
 文章快速检索 高级检索
 大地测量与地球动力学  2018, Vol. 38 Issue (10): 1033-1037,1067  DOI: 10.14075/j.jgg.2018.10.008

引用本文

SU Jinglan, ZHANG Hongping. Performance Analysis of RTK Algorithm for Single-Frequency Combination of GPS and BDS in Urban Environments[J]. Journal of Geodesy and Geodynamics, 2018, 38(10): 1033-1037,1067.

文章历史

1. 武汉大学卫星导航定位技术研究中心，武汉市珞喻路129号, 430079

1 BDS/GPS单频RTK算法模型

1.1 系统模型

 $\mathit{\boldsymbol{x}} = {\left( {\mathit{\boldsymbol{r}}_r^{\rm{T}},\mathit{\boldsymbol{v}}_r^{\rm{T}},c{\rm{d}}{{\mathit{\boldsymbol{\dot t}}}_r},\mathit{\boldsymbol{B}}_1^{\rm{T}}} \right)^{\rm{{\rm T}}}}$ (1)

 ${\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{r}}_r}}\\ {{\mathit{\boldsymbol{v}}_r}} \end{array}} \right)_t} = \underbrace {\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{3 \times 3}}}&{{\mathit{\boldsymbol{I}}_{3 \times 3}}\tau }\\ {{{\bf{0}}_{3 \times 3}}}&{{\mathit{\boldsymbol{I}}_{3 \times 3}}} \end{array}} \right)}_{{\mathit{\boldsymbol{F}}_{11}}}{\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{r}}_r}}\\ {{\mathit{\boldsymbol{v}}_r}} \end{array}} \right)_{t - 1}}$ (2)

GNSS单频接收机钟差随时间的变化较为频繁，数值差异也较为明显。当伪距、载波相位和多普勒观测值进行组合解算时，需要估计接收机的钟差变率，同时可以将估计得到的接收机钟差变率值用于辅助卫星钟跳、接收机钟跳的探测。考虑滤波方程的稳健性，文中将接收机钟差变率建模为一阶高斯-马尔科夫过程：

 $c{\rm{d}}{{\mathit{\boldsymbol{\ddot t}}}_r} = - \frac{1}{{{\tau _c}}}c{\rm{d}}{{\mathit{\boldsymbol{\dot t}}}_r} + \mathit{\boldsymbol{w}}$ (3)

 ${\left( {c{\rm{d}}{{\mathit{\boldsymbol{\dot t}}}_r}} \right)_t} \approx \underbrace {\exp \left( { - \frac{\tau }{{{\tau _c}}}} \right)}_{{\mathit{\boldsymbol{F}}_{22}}}{\left( {c{\rm{d}}{{\mathit{\boldsymbol{\dot t}}}_r}} \right)_{t - 1}}$ (4)

 ${\mathit{\boldsymbol{B}}_{1t}} = \underbrace {{\mathit{\boldsymbol{I}}_{m \times m}}}_{{\mathit{\boldsymbol{F}}_{33}}}{\mathit{\boldsymbol{B}}_{1t - 1}}$ (5)

 ${\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{r}}_r}}\\ {{\mathit{\boldsymbol{v}}_r}}\\ {c{\rm{d}}{{\mathit{\boldsymbol{\dot t}}}_r}}\\ {{\mathit{\boldsymbol{B}}_1}} \end{array}} \right)_t} = \left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{F}}_{11}}}&{}&{}\\ {}&{{\mathit{\boldsymbol{F}}_{22}}}&{}\\ {}&{}&{{\mathit{\boldsymbol{F}}_{33}}} \end{array}} \right){\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{r}}_r}}\\ {{\mathit{\boldsymbol{v}}_r}}\\ {c{\rm{d}}{{\mathit{\boldsymbol{\dot t}}}_r}}\\ {{\mathit{\boldsymbol{B}}_1}} \end{array}} \right)_{t - 1}}$ (6)
1.2 量测模型

BDS/GPS单频RTK利用流动站与基准站的伪距、载波相位和多普勒观测值组成双差观测值y =(Φ1T, P1T, D1T)T，以此更新滤波系统的量测方程h(x)=(hΦ, 1T, hP, 1T, hD, 1T)T。量测模型如式(7)所示：

 $\begin{array}{*{20}{c}} {\left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{rb}^i}\\ {\mathit{\boldsymbol{P}}_{rb}^i}\\ {\mathit{\boldsymbol{D}}_{rb}^i} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\rho }}_{rb}^i + {\mathit{\boldsymbol{\lambda }}_{1 * }}\mathit{\boldsymbol{B}}_{rb}^i}\\ {\mathit{\boldsymbol{\rho }}_{rb}^i}\\ {\mathit{\boldsymbol{p}}_r^i} \end{array}} \right) \approx {{\left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\rho }}_{rb}^i + {\mathit{\boldsymbol{\lambda }}_{1 * }}\mathit{\boldsymbol{B}}_{rb}^i}\\ {\mathit{\boldsymbol{\rho }}_{rb}^i}\\ {\mathit{\boldsymbol{p}}_r^i} \end{array}} \right)}_0} + }\\ {\left( {\begin{array}{*{20}{c}} { - \mathit{\boldsymbol{E}}}&{\bf{0}}&{\bf{0}}&{{\mathit{\boldsymbol{\lambda }}_{1 * }}{\mathit{\boldsymbol{I}}_{m \times m}}}\\ { - \mathit{\boldsymbol{E}}}&{\bf{0}}&{\bf{0}}&{\bf{0}}\\ {{\mathit{\boldsymbol{H}}_D}}&{ - \mathit{\boldsymbol{E}}}&{\bf{1}}&{\bf{0}} \end{array}} \right){\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{x}}} \end{array}$ (7)

 $\mathit{\boldsymbol{p}} = \mathit{\boldsymbol{\dot r}}_r^{\rm{s}} + c{\rm{d}}{{\mathit{\boldsymbol{\dot t}}}_r} - c{\rm{d}}\mathit{\boldsymbol{\dot T}}$ (8)
 $\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\dot r}}_r^{\rm{s}} = \mathit{\boldsymbol{e}}_r^{{\rm{sT}}}\left( {{\mathit{\boldsymbol{v}}^{\rm{s}}}\left( {{\mathit{\boldsymbol{t}}^{\rm{s}}}} \right) - {\mathit{\boldsymbol{v}}_r}} \right) + }\\ {\frac{{{\mathit{\boldsymbol{w}}_e}}}{c}\left( {\mathit{\boldsymbol{v}}_y^s{x_r} + {y^{\rm{s}}}{\mathit{\boldsymbol{v}}_{x,r}} - \mathit{\boldsymbol{v}}_x^s{y_r} - {x^{\rm{s}}}{\mathit{\boldsymbol{v}}_{y,r}}} \right)} \end{array}$ (9)

1.3 整数参数约束模型

 $\left( {\begin{array}{*{20}{c}} {\bf{0}}&{\bf{0}}&{\bf{0}}&{{\mathit{\boldsymbol{\lambda }}_{1 * }}}&{{\mathit{\boldsymbol{I}}_{m \times m}}} \end{array}} \right)\delta x = {\mathit{\boldsymbol{B}}_{t - k}} + {\mathit{\boldsymbol{v}}_b}$ (10)

1.4 随机模型

 $\sigma _{P/\mathit{\Phi }}^2 = \left\{ \begin{array}{l} {F_s}\left( {{a^2} + \frac{{{b^2}}}{{\sin e}}} \right) + C \times {10^{ - \frac{{C/N}}{{10}}}},{\rm{EL}} < {\rm{3}}{{\rm{0}}^ \circ }\\ {F_s}\left( {{a^2} + {b^2}} \right),{\rm{EL}} \ge {\rm{3}}{{\rm{0}}^ \circ } \end{array} \right.$ (11)

 ${\sigma _D} = \frac{{{\lambda ^ * }L}}{{2{\rm{ \mathsf{ π} }}T}}\sqrt {\frac{{4F{B_n}}}{{C/N}}\left( {1 + \frac{1}{{{T_D} \cdot C/N}}} \right)}$ (12)

2 BDS/GPS单频RTK模糊度固定

BDS/GPS单频RTK模糊度固定算法模型如下[10]

 $\mathit{\boldsymbol{y}} = N\left( {Aa + Bb,{P_{yy}}} \right),a \in {\mathbb{Z}^n},b \in {\mathbb{R}^p}$ (13)

 $\left( {\begin{array}{*{20}{c}} {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} }\\ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} } \end{array}} \right) \sim N\left( {\left[ {\begin{array}{*{20}{c}} a\\ b \end{array}} \right],\left[ {\begin{array}{*{20}{c}} {{P_{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} }}}&{{P_{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} }}}\\ {{P_{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} }}}&{{P_{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} }}} \end{array}} \right]} \right)$ (14)

$\vartheta$算子有多种形式，如LAMBDA/MLAMBDA、直接取整法(IR)、序贯取整法(IB)以及整数最小二乘求解法(ILS)等，本文模糊度固定采用MLAMBDA方法[11]

 ${\rm{ratio}} = \frac{{{{\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over a} }_2} - \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} } \right)}^{\rm{T}}}\mathit{\boldsymbol{P}}_a^{ - 1}\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over a} }_2} - \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} } \right)}}{{{{\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over a} }_1} - \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} } \right)}^{\rm{T}}}\mathit{\boldsymbol{P}}_a^{ - 1}\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over a} }_1} - \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} } \right)}} > {\mathit{\boldsymbol{R}}_{{\rm{thres}}}}$ (15)

 $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over b} = \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} - {\mathit{\boldsymbol{P}}_{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over b} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} }}\mathit{\boldsymbol{P}}_{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} }^{ - 1}\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a} - \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over a} } \right)$ (16)

 $\frac{{{{\left( {\mathit{\boldsymbol{y}} - h\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over a} }_2},{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over b} }_2}} \right)} \right)}^{\rm{T}}}\mathit{\boldsymbol{P}}_y^{ - 1}\left( {\mathit{\boldsymbol{y}} - h\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over a} }_2},{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over b} }_2}} \right)} \right)}}{{{{\left( {\mathit{\boldsymbol{y}} - h\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over a} }_1},{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over b} }_1}} \right)} \right)}^{\rm{T}}}\mathit{\boldsymbol{P}}_y^{ - 1}\left( {\mathit{\boldsymbol{y}} - h\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over a} }_2},{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over b} }_2}} \right)} \right)}} > {{\mathit{\boldsymbol{R'}}}_{{\rm{thres}}}}$ (17)

2.1 速度信息辅助模糊度固定

BDS/GPS模糊度的固定在很大程度上取决于浮点解的精度。在城市环境中伪距和载波相位观测值的噪声较大，一般浮点解的精度在1.0 m左右，然而采用多普勒观测值估计接收机运动速度，其精度在0.1 m/s左右[12]，当模糊度固定后，接收机运动速度精度也会有明显的提高。若充分利用速度信息以辅助模糊度参数的固定，则在一定程度上可以提高模糊度固定率。

 $\begin{array}{*{20}{c}} {{{\left( {{{\bar r}_r}} \right)}_t} = {{\left( {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over r} }_r}} \right)}_{t - M}} + }\\ {\sum\limits_{n = 1}^M {\left[ {v\left( {t - n + 1} \right) + v\left( {t - n} \right)} \right]\frac{{{T_n}}}{2}} } \end{array}$ (18)

 $\sigma _{{{\bar r}_r}}^2 = \sigma _{{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over r} }_{t - M}}}^2 + \sum\limits_{n = 1}^M {\left[ {\sigma _{v\left( {t - n + 1} \right)}^2 + \sigma _{v\left( {t - n} \right)}^2} \right]\frac{{T_n^2}}{4}}$ (19)

 $\left( {\begin{array}{*{20}{c}} 1&0&0&0 \end{array}} \right)\delta x = {\left( {{{\bar r}_r}} \right)_t} + v_r^ -$ (20)

2.2 模糊度继承策略

1) 维持一个双差模糊度列表用于保存固定成功的双差模糊度(B, nFix, dH, Age)，同时设置继承窗口大小K值、可信任继承阈值Q以及可继承时间限制Agemax(以s为单位)。其中，B表示模糊度; nFix表示的是模糊度固定次数，当模糊度固定时加1，反之减1;dH表示当前历元模糊度浮点解与列表中模糊度的差值; Age表示模糊度固定时间与当前时间的差值。K值和Q值的选取与测量噪声和速度有关，当测量噪声较大时，K取值较小，Q取值较大，本文采用Q=1.5，Agemax=3.0，并分析不同K值对单频RTK解算性能的影响。

2) 当MLAMBDA方法固定模糊度失败时，遍历列表查找满足条件(21)的双差模糊度：

 ${n_{{\rm{Fix}}}} \ge K,{d_H} \ge Q,{\rm{Age}} \le {\rm{Ag}}{{\rm{e}}_{\max }}$ (21)

3) 采用式(17)计算非模糊度参数的固定解，并用式(16)和式(18)检验模糊度是否继承成功。若模糊度继承成功，则更新列表。

3 实验分析

2017-07-13在武汉市区进行车载实验，设计路线包括林荫道、高楼建筑、高架桥、隔音墙等，其测试轨迹如图 1所示。测试中搭载了武汉迈普时空导航科技有限公司的战术级GNSS+INS组合导航系统POS310(陀螺零偏0.5°/h，加速度计零偏50 mGal)；配置了高精度Trimble BD982 GNSS OEM板卡，接收机采样率为1 Hz；基准站采用Trimble Net R9接收机，架设在武汉大学教学实验大楼楼顶，基线长度为11.6~16.3 km。由于城市环境下GNSS RTK/INS紧组合可以大大提高系统的可靠性和定位精度，并提供连续的高精度导航定位结果，故本文将NovAtel商业软件Inertial Explorer 8.4的GNSS RTK/INS双频紧组合结果作为轨迹的参考真值，并使用其中的单频数据进行解算，以此分析BDS/GPS单频RTK的定位精度。

 图 1 测试轨迹、典型场景以及GNSS接收机设备 Fig. 1 Vehicle trajectory, typical scenario and GNSS receiver equipment

 图 2 GPS、BDS和BDS+GPS单频RTK解算坐标偏差 Fig. 2 E, N, U bias of GPS, BDS and BDS+GPS single-frequency RTK solution

 图 3 GPS、BDS和BDS/GPS单频RTK解算结果误差统计 Fig. 3 Error statistics of GPS, BDS and BDS+GPS single-frequency RTK solution

 图 4 有/无速度辅助信息辅助时BDS、GPS和BDS+GPS单频RTK解算结果坐标偏差序列 Fig. 4 Bias of BDS, GPS and BDS+GPS single-frequency RTK solution with/without aid-informations of speeds

4 结语

 [1] 胡国荣, 欧吉坤. 改进的高动态GPS定位自适应卡尔曼滤波方法[J]. 测绘学报, 1999, 28(4): 290-294 (Hu Guorong, Ou Jikun. The Improved Method of Adaptive Kalman Filter for GPS High Kinematic Positioning[J]. Acta Geodaetica Cartographica Sinica, 1999, 28(4): 290-294 DOI:10.3321/j.issn:1001-1595.1999.04.003) (0) [2] 何海波, 李金龙, 郭海荣, 等. 北斗/GPS双系统单频RTK模糊度解算性能分析[J]. 测绘科学与工程, 2014, 34(1): 50-54 (He Haibo, Li Jinlong, Guo Hairong, et al. Performance Analysis of Single-Frequency RTK Ambigutity Resolution Using Dual System Beidou/GPS[J]. Geomatics Science and Engineering, 2014, 34(1): 50-54) (0) [3] Chui C K, Chen G. Kalman Filtering with Real-Time Applications[M]. Springer-Verlag, 1987 (0) [4] Zhao S, Cui X, Guan F, et al. A Kalman Filter-Based Short Baseline RTK Algorithm for Single-Frequency Combination of GPS and BDS[J]. Sensors, 2014, 14(8): 15415 DOI:10.3390/s140815415 (0) [5] Elsobeiey M E. Stochastic Analysis of Low-Cost Single-Frequency GPS Receivers[J]. Positioning, 2016, 7(3): 91-100 DOI:10.4236/pos.2016.73009 (0) [6] Odolinski R, Teunissen P J G. Low-Cost, High-Precision, Single-Frequency GPS-BDS RTK Positioning[J]. GPS Solutions, 2017, 21(1): 1-16 DOI:10.1007/s10291-016-0542-0 (0) [7] Stempfhuber W. 3D-RTK Capability of Single GNSS Receivers[C]. ISPRS-International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2013 (0) [8] Odolinski R, Teunissen P J G, Odijk D. Combined GPS, Beidou, Galileo, and QZSS Single-Epoch, Single-Frequency RTK Performance Analysis[C]. International Association of Geodesy Symposia, 2016 (0) [9] Odolinski R, Teunissen P J G, Odijk D. Combined BDS, Galileo, QZSS and GPS Single-Frequency RTK[J]. GPS Solutions, 2015, 19(1): 151-163 DOI:10.1007/s10291-014-0376-6 (0) [10] Takasu T, Yasuda A. Development of theLow-Cost RTK-GPS Receiver with an Open Source Program Package RTKLIB[C]. International Symposium on GPS/GNSS, South Korea, 2009 (0) [11] Chang X W, Yang X, Zhou T. MLAMBDA: A Modified LAMBDA Method for Integer Least-Squares Estimation[J]. Journal of Geodesy, 2005, 79(9): 552-565 DOI:10.1007/s00190-005-0004-x (0) [12] Kubo N. Advantage of Velocity Measurements on Instantaneous RTK Positioning[J]. GPS Solutions, 2009, 13(4): 271-280 DOI:10.1007/s10291-009-0120-9 (0)
Performance Analysis of RTK Algorithm for Single-Frequency Combination of GPS and BDS in Urban Environments
SU Jinglan1     ZHANG Hongping1
1. Research Center of GNSS, Wuhan University, 129 Luoyu Road, Wuhan 430079, China
Abstract: GPS/BDS dual-frequency RTK positioning is fast and reliable, and is most widely used in high-precision dynamic positioning. However, its receiver is expensive and difficult to be widely used in urban environments. It is easy to promote GNSS single-frequency RTK because of low-price GNSS receiver. In this paper, we introduce a single-frequency BDS+GPS RTK algorithm, taking account of velocity information to assist in ambiguity fix, the inheritance of ambiguity, and the fusion of Doppler observation. Through an actual on-board test in urban environment, the single-frequency RTK positioning performance of BDS, GPS and BDS+GPS are compared and analyzed.
Key words: urban environment; single-frequency RTK; integer-parameter constraint; velocity information assistance