兰勃脱等角圆锥投影反解不同算法的解析

丁士俊; 李鹏鹏; 邹进贵; 金银龙

doi:10.13203/j.whugis20200383

兰勃脱等角圆锥投影反解不同算法的解析

1.
武汉大学测绘学院，湖北武汉，430079
2.
武汉市测绘研究院，湖北武汉，430022
3.
武汉大学水利水电学院，湖北武汉，430072

基金项目:

国家自然科学基金 41871373

武汉大学“教育教学改革”建设引导专项 2022

详细信息

作者简介:
丁士俊，博士，教授，研究方向为现代测量数据处理理论、测量与地图投影理论与方法。shjding@sgg.whu.edu.cn

通讯作者:
李鹏鹏，博士。276149876@qq.com

中图分类号: P282
计量
- 文章访问数: 586
- HTML全文浏览量: 183
- PDF下载量: 52
出版历程
- 收稿日期: 2022-09-15
- 网络出版日期: 2022-09-19
- 发布日期: 2022-09-04

Analysis of Different Algorithms for Reverse Solution of Lambert Conic Conformal Projection

1.
School of Geodesy and Geomation, Wuhan University, Wuhan 430079, China
2.
Wuhan Geomatic Institute, Wuhan 430022, China
3.
School of Water Resources and Hydropower Engineering, Wuhan University, Wuhan 430072, China

Funds:

The National Natural Science Foundation of China 41871373

the Wuhan University Construction of Education and Teaching Reform 2022

More Information

Author Bio:
DING Shijun, PhD, professor, specializes in survey data processing, theories and methods of surveying and map projection.E-mail: shjding@sgg.whu.edu.cn

Corresponding author:
LI Pengpeng, PhD. E-mail: 276149876@qq.com

摘要

摘要: 在测量与地图制图中，等量纬度求解大地纬度是一种常见的投影反解计算，就该反解问题的几种不同算法进行研究，包括迭代法、等量纬差求解大地纬度的级数展开式及等量纬度求解大地纬度的直接算法。利用Mathematica对后两种算法的计算公式进行了详细推导，给出了其高阶系数展开式，同时对现有算法中存在的问题进行了解析。兰勃脱等角投影算例表明，所推导的公式其计算精度可达(1×10^-7)″~(1×10^-8)″，完全满足测量与地图投影高精度的要求。
- 等量纬度 /
- 大地纬度 /
- 投影反解 /
- 兰勃脱等角圆锥投影
Abstract:
Objectives In surveying and cartography, solving geodetic latitude is a common reverse calculation by conformal latitude.
Methods Several different algorithms of reverse solution are studied, including iterative method, series expansion of solving geodetic latitude by difference of conformal latitude and direct algorithm of inverse solution of geodetic latitude with conformal latitude. The calculation formulas of the latter two algorithms are derived by Mathematica. The higher-order coefficient expansion formula are given. At the same time, problems in existing algorithms are analyzed.
Results and Conclusions E xamples results of Lambert conic conformal projection show that the accuracy of the formula derived in this paper can reaches (1×10^-7)″-(1×10^-8)″. It fully meets the requirements of high precision in surveying and map projection.
- conformal latitude /
- geodetic latitude /
- reverse solution of projection /
- Lambert conic conformal projection

HTML全文

随着社会的进步和技术的发展，人们对于位置服务的需求日益增加。特别是自动驾驶、无人机和物联网等技术的日益进步，越来越多的人关注如何在复杂环境下实现精密导航定位。然而，任何单一的技术都不能够在所有环境下随时提供高精度、高可用性、高连续性和高可靠性的定位结果^[1]。因此，融合多种定位手段成为了导航定位领域研究的热点^[2-3]，全球卫星导航定位系统(global navigation satellite system，GNSS)、惯性导航系统(inertial navigation system，INS)和视觉的组合是其中一种具有潜力的组合方式。

GNSS具有全球、全天候、实时的导航定位授时能力，被广泛应用于交通等领域^[4]。Wübbena等^[5]于2005年提出PPP-RTK(precise point positioning real time kinematic)的概念，PPP-RTK急需解决的问题是模糊度的快速固定^[6]。文献[7-9]分别提出了未校正相位延迟(uncalibrated phase bias，UPD)、整数钟和解耦钟的固定方法，文献[10-11]基于实测数据对这3种方法进行了比较分析。近年来，PPP-RTK技术逐渐从无电离层组合向非差非组合方向发展。非差非组合通过约束电离层和对流层参数，在避免了无电离层组合将观测噪声放大了近3倍的缺点的同时，能进一步加速PPP的收敛^[12]。对非差非组合而言，通常认为大气增强和模糊度固定PPP技术即为PPP-RTK技术。PPP-RTK避免了RTK覆盖范围受限和PPP收敛时间长的缺陷，在车载导航等领域具有很广泛的应用前景。鉴于此，一些地区的官方机构已经推出了PPP-RTK服务，例如日本的准天顶星系统(quasi-zenith satellite system，QZSS)于2018年开始提供厘米级精度的服务^[13]。

尽管PPP-RTK具有众多优势，但是仍免不了受卫星信号被遮挡的影响，意味着在复杂环境下，PPP-RTK不能够提供高精度的定位服务。当信号受建筑物遮挡影响很严重时，可用卫星数少于4颗，此时PPP-RTK将不能提供定位结果。INS能在没有外部输入的情况下提供导航信息，但是误差随时间发散较快。INS在GNSS信号被遮挡时辅助提供定位结果，GNSS约束INS的误差发散，两者能形成良好的互补。Rabbou等^[14]和Gao等^[15]研究了PPP/INS紧组合，Gu等^[16]对大气增强PPP/INS紧组合进行了深入分析，但PPP-RTK/INS相关的研究相对较少，仍有待进一步深入研究。

PPP-RTK与INS进行组合在GNSS信号较长时间中断时，由于INS航位推算误差的发散，依旧产生较大的定位误差，因此需要进一步引入其他辅助的传感器来限制INS的误漂移误差。视觉里程计(visual odometry，VO)的误差发散较慢，可以对INS的漂移误差进行约束。但单目视觉里程计缺少尺度信息，因此可以通过INS和VO组合的方式，通过INS来恢复单目VO的尺度信息，即为视觉惯性里程计(visual inertial odometry，VIO)。VIO的融合方式主要分基于扩展卡尔曼滤波^[17-19]和非线性优化^[20-21]，前者结构简单，但线性化误差较大，后者精度更高，但对算力要求很高，计算负担大。多状态约束卡尔曼滤波(multistate constraint Kalman filter，MSCKF)^[17]是一种基于扩展卡尔曼滤波的VIO，能在大尺度环境下提供良好的状态估计。MSCKF将多帧影像对同一个特征点的观测转化为相机状态之间的约束，计算复杂度与特征点数仅成线性关系^[22]。

GNSS、INS和视觉的组合方式能更好地形成互补，GNSS为INS和视觉提供绝对位置信息，INS和视觉在GNSS信号被干扰时提供比较精确的航位推算信息。Kim等^[23]等使用六自由度的同步定位与测图(simultaneous localization and mapping，SLAM)辅助GNSS/INS在GNSS拒止位置环境下提供可靠的导航解，Won等^[24]提出了GNSS、INS和视觉的选择性组合方法，并且证明选择性组合方法要明显优于非组合的方式。以上这些组合中仅仅用到了GNSS的伪距定位，相位的潜力未能得到充分挖掘。Liu^[25]提出了GNSS/INS/立体视觉/地图匹配系统的紧组合模型；Li等^[26]实现了RTK、INS和单目相机的紧组合；Li等^[27]实现了城市复杂环境下PPP和立体视觉惯性导航系统(stereo visual INS，SVINS)的半紧组合，并进行了测试分析；方文涛^[28]实现了大气增强PPP/微机电系统(micro-electro-mechanical system，MEMS)惯导/视觉里程计融合定位研究。PPP和RTK与其他定位及技术的组合研究均有了一定的进展，但目前对于PPP-RTK的研究集中于改善PPP-RTK本身的性能，与其他定位技术组合的研究较少。PPP-RTK、INS和视觉的组合有望实现在城市环境下高精度连续定位，但其实际性能表现仍有待进一步研究。

对此，本文提出多系统非差非组合PPP-RTK与VIO的半紧组合模型，并且评估多系统PPP-RTK/VIO半紧组合在城市复杂环境下的定位表现，探究多系统PPP-RTK/VIO半紧组合在城市车载导航应用场景下的潜力。

1 PPP-RTK/VIO组合模型

PPP-RTK/VIO组合模型是由非差非组合PPP-RTK和单目VIO组成，基于非差非组合PPP-RTK和单目VIO模型推导，可得到组合模型。

1.1 PPP-RTK模型

非差非组合GNSS PPP观测模型通常表示如下^[29-30]：

\{\begin{array}{l} P_{r, f}^{s} = ρ_{r}^{s} + t_{r, s y s} + α_{r}^{s} T_{z} + \frac{40.3}{f^{2}} γ_{r}^{s} I_{r}^{s} - b^{s, f} + \\ b_{r, f} + ε_{P} \\ Φ_{r, f}^{s} = ρ_{r}^{s} + t_{r, s y s} + α_{r}^{s} T_{z} - \frac{40.3}{f^{2}} γ_{r}^{s} I_{r}^{s} + \\ λ_{f} N_{r, f}^{s} + ε_{Φ} \end{array}

式中， $P_{r, f}^{s}$ 、 $Φ_{r, f}^{s}$ 分别为伪距和相位观测值； $ρ_{r}^{s}$ 为卫星至接收机的几何距离； $t_{r, s y s}$ 为接收机钟差误差； $T_{z} 、 I_{r}^{s}$ 分别为天顶对流层误差和电离层误差； $α_{r}^{s}$ 、 $γ_{r}^{s}$ 分别为对应的投影函数； $b_{r, f} 、 b^{s, f}$ 分别为接收机端和卫星端硬件延迟误差； $N$ 、 $λ$ 分别为浮点模糊度和波长； $ε_{P} 、 ε_{Φ}$ 分别为伪距和相位观测噪声；角标 $f 、 s 、 r 和 s y s$ 分别表示对应的频率、卫星号、接收机 $和$ 系统，本文中主要涉及全球定位系统(global positioning system，GPS)和北斗卫星导航系统(BeiDou navigation satellite system，BDS)，分别简记为G和C。

对卫星端硬件延迟进行改正以及使用对流层的先验模型后，对式(1)进行线性化，则有：

\{\begin{matrix} Δ P_{r, f}^{s} = h_{r}^{s} δ x_{G N S S}^{e} + t_{r, s y s} + α_{r}^{s} δ T_{w} + (40.3 / f^{2}) γ_{r}^{s} I_{r}^{s} + b_{r, f} + ε_{P} \\ Δ Φ_{r, f}^{s} = h_{r}^{s} δ x_{G N S S}^{e} + t_{r, s y s} + α_{r}^{s} δ T_{w} - (40.3 / f^{2}) γ_{r}^{s} I_{r}^{s} + λ N_{r, f}^{s} + ε_{Φ} \end{matrix}

式中， $Δ P_{r, f}^{s}$ 、 $Δ Φ_{r, f}^{s}$ 分别为伪距和相位的验前残差(observed-minus-computed，OMC)； $h_{r}^{s}$ 为卫地距在接收机坐标系下的方向余弦向量； $δ x_{G N S S}^{e}$ 为接收机位置误差向量； $δ T_{w}$ 为天顶对流层湿延迟残差。对于电离层延迟的参数化，采用Shi等^[31]和Zhao等^[29]提出的DESIGN-5(deterministic plus stochastic ionosphere model for GNSS)模型，对电离层误差建模如下：

\{\begin{matrix} I_{r}^{s} = a_{0} + a_{1} d L + a_{2} d B + a_{3} d L^{2} + \\ a_{4} d B^{2} + r_{r}^{s} \\ {\tilde{I}}_{r}^{s} = a_{0} + a_{1} d L + a_{2} d B + a_{3} d L^{2} + \\ a_{4} d B^{2} + r_{r}^{s} + ε_{{\tilde{I}}_{r}^{s}} \end{matrix}

式中， $a_{i} (i = 0, 1, 2, 3, 4)$ 为描述电离层趋势项变化的系数，用于描述单站上空电离层延迟的空间分布，反映了电离层延迟的结构特性；随机参数 $r_{r}^{s}$ 用于表示电离层在时域上的随机特性；dL、dB为电离层穿刺点(ionospheric pierce point，IPP)与测站近似坐标的经纬度差值； ${\tilde{I}}_{r}^{s}$ 为电离层误差虚拟观测值，可以通过先验模型或者高精度电离层产品获取； $ε_{{\tilde{I}}_{r}^{s}}$ 是其对应的观测噪声。由此可将待估状态量 $x_{P P P - R T K}$ 表示为：

x_{P P P - R T K} = {(\begin{matrix} \begin{matrix} δ x_{G N S S}^{e} & t_{r} & δ T_{w} \end{matrix} & \begin{matrix} b_{r} & N_{r} & \begin{matrix} a_{r} & r_{r} \end{matrix} \end{matrix} \end{matrix})}^{T}

式中 $, t_{r} = t_{r, s y s}, s y s \in (G C)$ 为接收机对应GPS和BDS的钟差； $N_{r} = {(\begin{matrix} N_{r, 1} & N_{r, 2} \end{matrix})}^{T}$ 是两个频率上的浮点模糊度； $a_{r} = {(\begin{matrix} a_{0} & a_{1} & \begin{matrix} a_{2} & a_{3} & a_{4} \end{matrix} \end{matrix})}^{T}$ 和 $r_{r} = {(\begin{matrix} r_{r}^{1} r_{r}^{2} & \dots & r_{r}^{j} \end{matrix})}^{T}$ 是DESIGN模型的参数。

对于PPP-RTK中模糊度整数解，浮点模糊度 $N_{r, f}^{s}$ 进一步建模为：

N_{r, f}^{s} = n - d_{r} + d^{s}

式中，n代表整数模糊度； $d_{r}$ 、 $d^{s}$ 分别代表接收机和卫星端UPD(uncalibrated phase delays)。由式(2)得到浮点模糊度之后，通过服务端提供的UPD产品来恢复模糊度的整数特性。进一步采用LAMBDA(least-squares ambiguity decorrelation adjustment)^[32]方法来搜索最优的模糊度固定值。将模糊度固定值作为约束条件，得到模糊度固定的PPP定位结果^[33-34]。

1.2 INS模型

选用e系(地心地固系)作为机械编排的参考系，惯性系统的误差状态建模为：

x_{I N S} = {(\begin{matrix} \begin{matrix} δ x_{I N S}^{e} & δ v_{I N S}^{e} & ϕ \end{matrix} & \begin{matrix} B_{a} & B_{g} \end{matrix} & \begin{matrix} S_{a} & S_{g} \end{matrix} \end{matrix})}^{T}

式中， $δ x_{I N S}^{e}$ 表示地心地固坐标系下INS中心的位置改正向量； $δ v_{I N S}^{e}$ 表示地心地固坐标系下INS中心的速度改正向量； $ϕ$ 表示姿态误差； $B_{a}$ 、 $B_{g}$ 分别表示加速度计和陀螺仪的零偏误差； $S_{a}$ 、 $S_{g}$ 分别表示加速度计和陀螺仪的比例因子误差。对于上述的位置、速度、姿态误差，姿态误差模型采用 $ϕ$ 角模型，其连续型动态模型通常表达如下^[35]：

$$\begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {{\rm{\delta }}\dot x_{{\rm{INS}}}^e{\rm{\delta }}\dot v_{{\rm{INS}}}^e\dot \phi } \right)^{\rm{T}}} = \\ \left( {\begin{array}{*{20}{c}} {{\rm{\delta }}v_{{\rm{INS}}}^e}\\ { - {\rm{2}}\omega _{ie}^e \times {\rm{\delta }}v_{{\rm{INS}}}^e + R_b^e{f^b} \times \phi + R_b^e{\rm{\delta }}{f^b} + {\rm{\delta }}{g^e}}\\ { - \omega _{ie}^e \times \phi - R_b^e\delta \omega _{ib}^b} \end{array}} \right) \end{array}$$

式中， $ω_{i e}^{e}$ 为地球自转的角速度矢量； $R_{b}^{e}$ 为载体坐标系转换到地心地固系的方向余弦矩阵； $f^{b}$ 为经误差补偿后加速度计输出的三维比力矢量； $δ f^{b}$ 和 $δ ω_{i b}^{b}$ 是惯导的器件误差，考虑利用惯导零偏误差和比例因子误差来建模惯性器件误差方程：

\{\begin{matrix} δ ω_{i b}^{b} = B_{g} + d i a g (ω_{i b}^{b}) S_{g} + w_{g} \\ δ f^{b} = B_{a} + d i a g (f_{}^{b}) S_{a} + w_{a} \end{matrix}

式中 $, w_{g}$ 、 $w_{a}$ 分别是陀螺仪和加速度计的量测噪声。而对陀螺仪和加速度计零偏、比例因子误差建模为一阶高斯马尔可夫过程：

\{\begin{matrix} (\begin{matrix} {\dot{B}}_{a} \\ {\dot{B}}_{g} \end{matrix}) = (\begin{matrix} \frac{- 1}{τ_{b_{a}}} B_{a} \\ \frac{- 1}{τ_{b_{g}}} B_{g} \end{matrix}) + (\begin{matrix} w_{b_{a}} \\ w_{b_{g}} \end{matrix}) \\ (\begin{matrix} {\dot{S}}_{a} \\ {\dot{S}}_{g} \end{matrix}) = (\begin{matrix} \frac{- 1}{τ_{s_{a}}} S_{a} \\ \frac{- 1}{τ_{s_{g}}} S_{g} \end{matrix}) + (\begin{matrix} w_{s_{a}} \\ w_{s_{g}} \end{matrix}) \end{matrix}

式中， $w_{b_{a}} 、 w_{b_{g}} 、 w_{s_{a}}$ 、 $w_{s_{g}}$ 分别为加速度计和陀螺仪零偏误差、比例因子误差一阶高斯马尔可夫过程的驱动噪声； $τ_{b_{a}} 、 τ_{b_{g}}$ 、 $τ_{s_{a}} 、 τ_{s_{g}}$ 为对应的相关时间。

1.3 VIO模型

基于INS误差状态量，构建视觉惯导误差状态量。对第i时刻的相机状态而言，建立其误差状态为 $(δ p_{c_{i}}^{e} ϕ_{c_{i}})$ ，因此对第k时刻的基于MSCKF的VIO误差状态量建模如下：

x_{i, c} = {(\begin{matrix} x_{I N S} & \begin{matrix} \begin{matrix} | & δ p_{c_{j}}^{e} & ϕ_{c_{j}} \end{matrix} & \begin{matrix} \dots & δ p_{c_{k}}^{e} & ϕ_{c_{k}} \end{matrix} \end{matrix} \end{matrix})}^{T}

式中， $δ p_{c_{i}}^{e}$ 、 $ϕ_{c_{i}}$ 分别为i时刻(i = j，j+1，j+2… k)相机的位置误差向量和失准角向量，当有新时刻的视觉数据引入时，就对上述状态进行状态增广。

基于不在当前帧所观测到的静止特征点，可以构建多状态约束滤波的量测模型。通过数学变换实现将前n个时刻的相机状态与相应的静止特征点之间的约束关系转换为前n个时刻相机状态之间的约束关系。假定k+1时刻静止特征点 $P_{k}$ 不被观测到，但在 $P_{k}$ 被观测到的任意i时刻均有：

p_{P_{k}, i}^{c_{i}} = R_{e}^{c_{i}} (p_{P_{k}}^{e} - p_{c_{i}}^{e})

式中， $p_{P_{k}, i}^{c_{i}}$ 为 $P_{k}$ 点在第i时刻相机坐标系下的位置向量； $R_{e}^{c_{i}}$ 和 $p_{c_{i}}^{e}$ 分别表示第i时刻相机相对于e系的方向余弦矩阵和位置向量； $p_{P_{k}}^{e}$ 为 $P_{k}$ 点在e系下的位置向量，可通过多张影像三角测量获得。

构建相机量测误差状态量如下：

z_{P_{k}, i} = (\begin{matrix} u_{P_{k}, i}^{0} - {\tilde{u}}_{P_{k}, i} \\ v_{P_{k}, i}^{0} - {\tilde{v}}_{P_{k}, i} \end{matrix}), ε_{P_{k}, i} = (\begin{matrix} ε_{{\tilde{u}}_{P_{k}, i}} \\ ε_{{\tilde{v}}_{P_{k}, i}} \end{matrix})

式中， ${[u_{P_{k}, i}^{0}, v_{P_{k}, i}^{0}]}^{T}$ 为特征点 $P_{k}$ 重投影得到的像素坐标估计值； ${[{\tilde{u}}_{P_{k}, i}, {\tilde{v}}_{P_{k}, i}]}^{T}$ 为特征点 $P_{k}$ 在第i时刻的像素观测值； $ε_{P_{k}, i}$ 为对应的观测噪声。针对量测方程的系数阵，采用链式法则进行求解可得：

$$\begin{array}{l} {z_{{P_k},i}} = {\rm{\delta }}{\left( {{u_{{P_k},i}},{v_{{P_k},i}}} \right)^{\rm{T}}} = \frac{{\partial {{\left( {{u_{{P_k},i}},{v_{{P_k},i}}} \right)}^{\rm{T}}}}}{{\partial p_{{P_k},i}^{{c_i}}}}{\rm{\delta }}p_{{P_k},i}^{{c_i}} = \\ \;\;\;\;\;\;\;\;\;\;{H_{x,i}}{x_{i,c}} + {H_{f,i}}{\rm{\delta }}p_{{P_k}}^e + {\varepsilon _{{P_k},i}} \end{array}$$

式中，

H_{x, i} = [\begin{matrix} 0_{2 \times 21} & 0_{2 \times 6} & \dots & - J R_{e}^{c_{i}} & - J R_{e}^{c_{i}} [(p_{P_{k}}^{e} - p_{c_{i}}^{e}) \times] & \dots \end{matrix}]

H_{f, i} = J R_{e}^{c_{i}}

J = \frac{\partial {(u_{P_{k}, i}, v_{P_{k}, i})}^{T}}{\partial p_{P_{k}, i}^{c_{i}}} = (\begin{matrix} \frac{f_{x}}{z_{P_{k}, i}^{c_{i}}} & 0 & - \frac{f_{x} x_{P_{k}, i}^{c_{i}}}{{z_{P_{k}, i}^{c_{i}}}^{2}} \\ 0 & \frac{f_{y}}{z_{P_{k}, i}^{c_{i}}} & - \frac{f_{y} y_{P_{k}, i}^{c_{i}}}{{z_{P_{k}, i}^{c_{i}}}^{2}} \end{matrix})

$其中, (f_{x}, f_{y})$ 为相机焦距； $(\begin{matrix} x_{P_{k}, i}^{c_{i}} & y_{P_{k}, i}^{c_{i}} & z_{P_{k}, i}^{c_{i}} \end{matrix})$ 是特征点在相机坐标系下的位置；×表示由向量构建反对称矩阵。

由于 $δ p_{P_{k}}^{e}$ 不是待估状态量，并且对 $H_{f, i}$ 的量均为已知，这些量可由惯性测量单元(inertial measurement unit，IMU)机械编排获得，因此可以针对 $H_{f, i}$ 求解其左零空间矩阵A，使得：

A^{T} H_{f, i} = 0

对式(13)两边同时乘以 $A^{T}$ 矩阵，可得：

z_{a, i} = A^{T} z_{P_{k}, i} =

A^{T} H_{x, i} x_{i, c} + A^{T} ε_{P_{k}, i} = H_{a, x, i} x_{i, c} + ε_{a, i}

式中， $H_{a, x, i} = A^{T} H_{x, i}$ ； $ε_{a, i} = A^{T} ε_{P_{k}, i}$ 。

1.4 PPP-RTK/VIO半紧组合模型

在实际硬件平台上，惯导的位置 $x_{I N S}^{e}$ 与GNSS的天线相位中心位置 $x_{G N S S}^{e}$ 不重合，两者之间的空间关系表示为：

x_{G N S S}^{e} = x_{I N S}^{e} + C_{b}^{e} l^{b}

式中， $l^{b}$ 为杆臂信息，描述的是GNSS天线相位中心到IMU中心的向量在IMU坐标系下的投影。对式(16)进行微分可得到：

δ x_{G N S S}^{e} = δ x_{I N S}^{e} + C_{b}^{e} l^{b} \times ϕ

将PPP-RTK获得的位置信息记为 ${\tilde{x}}_{G N S S}^{e}$ ，VIO递推得到的位置记为 ${\hat{x}}_{G N S S}^{e}$ ，则可得到PPP-RTK/VIO的半紧组合观测方程：

\{\begin{array}{l} {\hat{x}}_{G N S S}^{e} - {\tilde{x}}_{G N S S}^{e} = δ x_{I N S}^{e} + C_{b}^{e} l^{b} \times ϕ + ε_{P} \\ z_{a, i} = H_{a, x, i} x_{i, c} + ε_{a, i} \end{array}

1.5 模型结构

PPP-RTK/VIO半紧组合算法结构如图 1所示，GNSS接收到观测数据后，结合大气增强数据进行PPP解算，得到浮点解后，加入UPD产品进行模糊度固定得到PPP-RTK定位结果。IMU进行机械编排推算下一时刻的位置，当接收到相机数据后，进行特征点提取和跟踪，相机和IMU数据同步则进行相机状态增广，并进行融合获得定位结果，当VIO与GNSS数据同步后，加入PPP-RTK定位结果到VIO滤波器中，得到PPP-RTK/VIO半紧组合输出的导航信息。

图 1 PPP-RTK/VIO半紧组合算法结构

Figure 1. Structure of Semi-Tightly Coupled PPP-RTK/VIO

下载: 全尺寸图片幻灯片

2 实验与结果分析

为了评估多系统PPP-RTK/VIO半紧组合的定位精度和可靠性，开发了FUSING(Fusing in GNSS)软件平台。目前，该平台实现了多系统轨道钟差估计，电离层、对流层建模、多频多系统偏差解算以及多源融合定位数据处理等^[36]。本文在FUSING平台的基础上，进一步实现了多系统PPP-RTK/VIO半紧组合数据处理，并利用车载实验进行了验证。

2.1 实验说明

测试数据所采用的车载平台如图 2所示。MEMS级惯导采集的数据适用于测试评估组合系统的表现，导航级惯导适用于计算参考解。两种IMU采样频率都是200 Hz，设备参数如表 1所示。灰度相机Basler acA640-90gm以659×494的分辨率、10 Hz的频率采集图像数据。UBLOX-M8T用于时间同步，产生PPS(pulse per second)信号触发相机曝光。GNSS数据由Trimble Net R9接收机以1 Hz的频率采集。IMU中心和GNSS天线相位中心的杆臂补偿在实验前手动测量，IMU和相机的外参由联合标定工具Kalibr( https://github.com/ethz-asl/kalibr/)标定。

图 2 GNSS/INS/视觉数据采集平台

Figure 2. GNSS/INS/Vision Data Collection Platform

下载: 全尺寸图片幻灯片

表 1 惯导设备参数

Table 1. Performance Specification of the IMU Sensors

惯导设备	随机游走		零偏稳定性
惯导设备	速度/( $m \cdot s$ ^－1 $\cdot$ $\sqrt[]{h}$ ^－1)	角度/( $° \cdot \sqrt[]{h}$ ^－1)	陀螺/( $° \cdot h$ ^－1)	加速度计/mGal
MEMS级惯导	0.03	0.17	8	200
导航级惯导	0.003	0.03	0.027	15

下载: 导出CSV

| 显示表格

测试数据为2020-01-02于武汉大学校内采集，采集数据时车速大约为5 m/s。图 3所示为测试轨迹和真实场景，可以看到测试环境中遮挡物较多，道路狭窄，两侧道路上都有较多的树木，还有一些建筑物，易对GNSS信号造成干扰。

图 3 测试轨迹和典型场景

Figure 3. Test Trajectory and Typical Scenarios

下载: 全尺寸图片幻灯片

如图 4所示为可用卫星数和截止角为10°时的位置精度强弱度(position dilution of precision，PDOP)值，GPS和GPS+BDS的平均可用卫星数分别为5.21颗和10.72颗。从图 4可以看出可用卫星数波动较大，说明GNSS信号受干扰比较严重；GPS和GPS+BDS的PDOP值分别为20.06和4.45，可见，卫星可用数下降时，PDOP值显著增大。

图 4 测试环境下的可用卫星数和PDOP值

Figure 4. Satellite Number of GPS/BDS and PDOP in Test Environment

下载: 全尺寸图片幻灯片

考虑到解算大气增强和UPD产品需要基站数据，收集了武汉周边的基站数据，图 5所示为用于解算大气增强和UPD产品的基站分布图，一共有7个基站，站间平均距离为40 km。图 5中红色轨迹为车载测试的轨迹。GNSS数据处理策略如表 2所示，定位结果以RMS(root mean square)的方式评估，用RTK/INS双向平滑松组合作为参考解。参考解采用商业软件GINS计算的结果，其标称精度可达到平面2 cm，高程3 cm。

图 5 基站分布图

Figure 5. Distribution of Reference Stations

下载: 全尺寸图片幻灯片

表 2 GNSS数据处理策略

Table 2. GNSS Data Processing Strategy

相关参数	PPP	PPP-RTK
模糊度	浮点解	固定解
对流层	GPT2w模型和VMF1_HT^[37], 残差建模为随机游走	区域模型
电离层	DESIGN-5模型^{[29, 31]}，GIM作为先验约束	区域模型
系统	G	G/G+C
观测值类型	非差非组合
观测值频点	GPS系统为C1C、C2W，北斗系统为C2I、C7I
PCO/PCV	采用igs14.atx改正
截止角/(°)	10
固体潮	IERS 2010
星历	GFZ精密星历产品
接收机钟差	建模为白噪声
未矫正码延迟	接收机端建模为随机游走，卫星端采用IGS产品直接改正
伪距观测值权重/m	0.3
相位观测值权重/m	0.003

下载: 导出CSV

| 显示表格

2.2 结果分析

所有解算的结果均以1 Hz的频率输出。G、C分别代表GPS系统和北斗系统，LC代表松组合。从图 6(a)可以看出，G-PPP-RTK明显优于G-PPP的定位效果，加入北斗系统后，定位误差序列的波动减小，定位精度更高，但是仍然是跳跃很明显，定位不够稳定。结合图 4卫星可用数可以看出，单GPS有很多历元可用卫星数小于4颗，此时无法解算出定位结果，定位误差序列图中被GC-PPP-RTK的定位误差序列掩盖而难以直观看出。但从东向周内秒355 250 s处(图 6(a)中黑框所示)仍旧可以明显看出G-PPP定位中断，而北斗的引入改善了定位连续性。加入惯导后(图 6(b))，GC-PPP-RTK因受信号干扰而定位精度不稳定的情况得到明显改善，水平和高程方向均很明显。相比于加入惯导，直接引入VIO，对GC-PPP-RTK的精度提升效果更加明显，北、东、地3个方向的定位误差序列都趋于比较稳定的状态，由此可以看出，GC-PPP-RTK/VIO的组合方式明显优于其余4种定位解算策略。

图 6 不同解算策略的定位误差时序图

Figure 6. Position Difference Series of Different Solutions

下载: 全尺寸图片幻灯片

为了进一步验证结论，表 3中列出了各种解算策略的统计精度。其中可用性定义为输出定位结果的历元数与该时段内的总历元数的比值，以此来评估定位的连续性。从表 3可以看出，北斗系统的引入大大提高了定位的连续性。定位精度方面，由于定位连续性显著提高，GC-PPP-RTK在信号被干扰的情况下输出定位结果的历元数增加，导致误差较大的定位结果占比增加，因此高程方向上的统计结果要比G-PPP和G-PPP-RTK差。为了体现不同组合方式的定位效果，表 4进一步以GC-PPP-RTK为参考，给出了GC-PPP-RTK/INS松组合和GC-PPP-RTK/VIO半紧组合对定位精度的提升效果。统计结果显示，引入惯导、VIO后，水平和高程相对于GC-PPP-RTK，均带来了20%~50%的精度提升，GC-PPP-RTK/VIO的定位效果最好，定位精度水平达到0.58 m，高程达到1.12 m，基本满足车道级的定位需求。

表 3 PPP和PPP-RTK的定位误差(RMS)统计结果

Table 3. RMS of the Position Difference of PPP and PPP-RTK

解算策略	水平/m	精度提升/%	高程/m	精度提升/%	可用性/%
G-PPP	0.97		2.21		82.2
G-PPP-RTK	0.91	6.2	1.46	33.9	82.2
GC-PPP-RTK	0.85	12.4	2.29	－3.6	99.4

下载: 导出CSV

| 显示表格

表 4 不同组合方式的定位误差(RMS)统计结果

Table 4. RMS of the Position Difference of Different Integration Models

解算策略	水平/m	精度提升/%	高程/m	精度提升/%	可用性/%
GC-PPP-RTK	0.85		2.29		99.4
GC-PPP-RTK/INS LC	0.62	27.1	1.64	28.4	100
GC-PPP-RTK/VIO	0.58	31.8	1.12	51.1	100

下载: 导出CSV

| 显示表格

考虑到在车载导航中速度和姿态也是重要的导航信息，本文也对组合系统的测速和姿态估计能力进行了评估。以GC-PPP-RTK/INS LC作为参考解，从图 7可以看出，GC-PPP-RTK/VIO的测速精度明显更好，精度更高，也更加稳定。由于GC-PPP-RTK高程方向定位误差较大且波动非常明显，一定程度上导致GC-PPP-RTK/INS LC高程方向的测速精度较低且波动很大，明显低于平面测速精度。GC-PPP-RTK/VIO组合方式对高程方向的定位约束较好，定位精度较高，高程方向的测速精度也明显提高。从表 5可以看出，GC-PPP-RTK/VIO在北向、东向和地向的测速精度分别达到0.04 m/s、0.04 m/s和0.02 m/s，相比于GC-PPP-RTK/INS LC，在地向的速度估计提升了83.3%，效果非常显著，北向和东向也分别提升了33.3%和20.0%。

图 7 不同解算策略的测速误差时序图

Figure 7. Velocity Estimation Difference Series of Different Solutions

下载: 全尺寸图片幻灯片

表 5 不同解算策略的测速误差(RMS)统计结果

Table 5. RMS of the Velocity Estimation Difference of Different Solutions

解算策略	北向/(m·s^－1)	东向/(m·s^－1)	地向/(m·s^－1)
GC-PPP-RTK/INS LC	0.06	0.05	0.12
GC-PPP-RTK/VIO	0.04	0.04	0.02
精度提升/%	33.3	20.0	83.3

下载: 导出CSV

| 显示表格

由于车载导航中，横滚角和俯仰角都是小角度且比较稳定，所以两种解算策略对横滚角和俯仰角的估计精度没有明显差别，这从图 8的误差序列和表 6的统计结果都可以看出来。而对于航向角，两种解算策略的估计精度则差别明显，GC-PPP-RTK/VIO对航向角的估计精度明显优于GC-PPP-RTK/INS LC，提升达到73.0%。

图 8 不同解算策略的姿态估计误差时序图

Figure 8. Pose Estimation Difference Series of Different Solutions

下载: 全尺寸图片幻灯片

表 6 不同解算策略的姿态估计误差(RMS)统计结果

Table 6. RMS of the Pose Estimation Difference of Different Solutions

解算策略	横滚角/(°)	俯仰角/(°)	航向角/(°)
GC-PPP-RTK/INS LC	0.10	0.06	0.63
GC-PPP-RTK/VIO	0.10	0.06	0.17
精度提升/%	0	0	73.0

下载: 导出CSV

| 显示表格

3 结语

测试结果表明，为了改善城市环境下车载导航的定位定姿和测速表现，本文提出了多系统PPP-RTK/VIO半紧组合算法模型，并利用城市环境车载数据进行了测试，测试环境选用武汉大学校园，为典型的GNSS信号受干扰比较严重的场景。另外收集了武汉周边的基站数据用于解算大气增强和UPD产品。

城市复杂环境下，PPP-RTK定位精度相比于PPP有明显改善，但定位精度仍然有待提高。PPP-RTK与VIO组合能带来明显的精度提高，相比于GC-PPP-RTK，平面和高程定位精度分别提升了31.8%和51.1%，达到了平面0.58 m、高程1.12 m的定位精度水平，基本满足车道级定位的需求。在测速和姿态估计方面，GC-PPP-RTK/VIO明显优于GC-PPP-RTK/INS LC，测速精度北向、东向和地向分别达到了0.04 m/s、0.04 m/s和0.02 m/s，横滚角、俯仰角和航向角的估计精度分别达到了0.10°、0.06°和0.17°。

本文实现了多系统PPP-RTK/VIO半紧组合算法，PPP-RTK/VIO紧组合算法模型仍有待进一步研究。PPP-RTK/VIO紧组合算法在使用一个滤波器的情况下，如何维持滤波器的稳定性，VIO的加入能否给PPP-RTK的模糊度固定带来明显的改善效果，以及其定位定姿和测速潜力等，都有待进一步研究和挖掘。

表 1 兰勃脱等角割圆锥投影正解

Table 1 Forward Calculation for Lambert Conic Conformal Projection with Two Standard Parallels

点号	B	L	X/m	Y/m
1	23°30'25.369 43″	46°50'47.284 55″	2 946 710.860 07	1 188 342.791 39
2	28°00'45″	45°30'25″	3 444 391.831 74	1 049 914.561 05

下载: 导出CSV

表 2 兰勃脱等角割圆锥投影反解

Table 2 Reverse Calculation for Lambert Conic Conformal Projection with Two Standard Parallels

点号	迭代法（ε=5×10^-14 rad）		等量纬差法	直接法
1	B	23°30'25.369 429 99″	23°30'25.369 429 99″	23°30'25.369 429 76″
	L	46°50'47.284 550 09″	46°50'47.284 550 09″	46°50'47.284 550 09″
	\|dB\|	（1×10^-8）″	（1×10^-8）″	（2.4×10^-7）″
2	B	28°00'45.000 000 034″	28°00'45.000 000 035″	28°00'44.999 999 945″
	L	45°30'25.000 000 121″	45°30'25.000 000 122″	45°30'25.000 000 122″
	\|dB\|	（3.4×10^-8）″	（3.5×10^-8）″	（4.5×10^-8）″

下载: 导出CSV

表 3 等量纬差反解大地纬度的比较

Table 3 Comparison of Reverse Solutions to Geodetic Latitude byConformal Latitude Difference

计算方法	1号点		2号点		系数项
计算方法	B	\|dB\|	B	\|dB\|	系数项
文献[12]算法（式（37））	23°30'25.369 433 56″	（3.6×10^-6）″	28°00'45.368 167 85″	0.681 68″	^*t₁~t₅
文献[12]算法（式（37））	23°30'25.369 330 48″	（9.9×10^-5）″	28°00'44.505 910 06″	0.494 08″	^*t₁~t₃
本文算法	23°30'25.369 429 99″	（1.0×10^-8）″	28°00'45.000 000 035″	（3.5×10^-8）″	t₁~t₅
本文算法	23°30'25.369 429 32″	（6.8×10^-7）″	28°00'45.999 998 692″	（1.3×10^-6）″	t₁~t₃
注：带^*的为文献[12]给出的系数

下载: 导出CSV

参考文献(21)

[1]	杨启和. 地图投影变换原理与方法[M]. 北京: 解放军出版社, 1989 Yang Qihe. The Principle and Method of Map Projection Transformation[M]. Beijing: The People's Liberation Army Press, 1989
[2]	杨启和. 测量和地图学中应用的六种纬度及其变换关系式[J]. 测绘科技通讯, 1995(3): 14-19 https://www.cnki.com.cn/Article/CJFDTOTAL-CHKJ199503004.htm Yang Qihe. Six Latitudes and Their Transformation Relations Used in Surveying and Cartography[J]. Geomatics Technology and Equipment, 1995(3): 14-19 https://www.cnki.com.cn/Article/CJFDTOTAL-CHKJ199503004.htm
[3]	李厚朴, 边少锋, 陈良友. 等面积纬度函数和等量纬度变换的直接解算公式[J]. 武汉大学学报·信息科学版, 2011, 36(7): 843-846 http://ch.whu.edu.cn/article/id/588 Li Houpu, Bian Shaofeng, Chen Liangyou. The Direct Calculating Formulae for Transformations Between Authalic Latitude Function and Isometric Latitude[J]. Geomatics and Information Science of Wuhan University, 2011, 36(7): 843-846 http://ch.whu.edu.cn/article/id/588
[4]	李厚朴, 边少锋, 刘敏. 地图投影中三种纬度间变换直接展开式[J]. 武汉大学学报·信息科学版, 2013, 38(2): 217-220 http://ch.whu.edu.cn/article/id/6108 Li Houpu, Bian Shaofeng, Liu Min. Direct Expansions of Transformations Between Three Kinds of Latitudes Used in Map Projection[J]. Geomatics and Information Science of Wuhan University, 2013, 38(2): 217-220 http://ch.whu.edu.cn/article/id/6108
[5]	Snyder J P. Map Projections Used by the U S Geological Survey[M]. Washington, USA: United States Government Printing Office, 1982
[6]	Snyder J P. Map Projections―A Working Manual[M]. Washington, USA: United States Government Printing Office, 1987
[7]	Grafarend E W, You R J, Syffus R. Map Projections[M]. Berlin, Germany: Springer, 2014
[8]	杨启和, 杨晓梅. 测量和地图学中应用的三种纬度函数及其反解变换的线性插值方法[J]. 测绘学报, 1997, 26(1): 92-94 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB701.016.htm Yang Qihe, Yang Xiaomei. Linear Interpolation Method of Three Latitude Functions and Their Inverse Solutions Used in Surveying Cartography[J]. Acta Geodaetica et Cartographic Sinica, 1997, 26(1): 92-94 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB701.016.htm
[9]	郑彤, 边少锋. 子午线弧长反问题新解[J]. 武汉大学学报·信息科学版, 2007, 32(3): 255-258 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200703016.htm Zheng Tong, Bian Shaofeng. New Solutions on Inverse Problem of Meridian Arc[J]. Geomatics and Information Science of Wuhan University, 2007, 32(3): 255-258 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200703016.htm
[10]	Craig R. Auxiliary Latitude Formulas: Finding the Coefficients Numerically and Symbolically[C]//Wolfram Technology Conference, Champaign, USA, 2006
[11]	李国藻, 杨启和, 胡定荃. 地图投影[M]. 北京: 解放军出版社, 1993 Li Guozao, Yang Qihe, Hu Dingquan. Map Projections [M]. Beijing: The People's Liberation Army Press, 1993
[12]	孔祥元, 郭际明, 刘宗泉. 大地测量学基础[M]. 武汉: 武汉大学出版社, 2010 Kong Xiangyuan, Guo Jiming, Liu Zongquan. Foundation of Geodesy[M]. Wuhan: Wuhan University Press, 2010
[13]	International Association of Oil & Gas Producers. Coordinate Conversions and Transformations Including Formulas[EB/OL]. [2018-01-25]. http//www. epsg. org/guides/index. html
[14]	边少锋, 纪兵. 等距离纬度等量纬度和等面积纬度展开式[J]. 测绘学报, 2007, 36(2): 218-223 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB200702017.htm Bian Shaofeng, Ji Bing. The Expansions of Rectifying Latitude, Conformal Latitude and Authalic Latitude[J]. Acta Geodaetica et Cartographica Sinica, 2007, 36(2): 218-223 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB200702017.htm
[15]	王瑞, 李厚朴. 辅助纬度反解公式的Lagrange级数法推演[J]. 海洋测绘, 2008, 28(3): 18-23 https://www.cnki.com.cn/Article/CJFDTOTAL-HYCH200803007.htm Wang Rui, Li Houpu. The Derivation of the Inverse Expansions for Auxiliary Latitudes by Lagrange Series Method[J]. Hydrographic Surveying and Charting, 2008, 28(3): 18-23 https://www.cnki.com.cn/Article/CJFDTOTAL-HYCH200803007.htm
[16]	李厚朴, 边少锋. 辅助纬度反解公式的Hermite插值法新解[J]. 武汉大学学报·信息科学版, 2008, 33(6): 623-626 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200806020.htm Li Houpu, Bian Shaofeng. Derivation of Inverse Expansions for Auxiliary Latitudes by Hermite Interpolation Method[J]. Geomatics and Information Science of Wuhan University, 2008, 33(6): 623-626 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200806020.htm
[17]	过家春, 李厚朴, 庄云玲, 等. 依不同纬度变量的子午线弧长正反解公式的级数展开[J]. 测绘学报, 2016, 45(5): 560-565 Guo Jiachun, Li Houpu, Zhuang Yunling, et al. Series Expansion for Direct and Inverse Solutions of Meridian in Terms of Different Latitude Variables[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(5): 560-565
[18]	陈成, 边少锋, 刘强. 纬度正反解问题的傅里叶级数解算方法[J]. 海洋测绘, 2017, 37(4): 19-23 Chen Cheng, Bian Shaofeng, Liu Qiang. A Fourier Series Method for Forward and Inverse Solutions of Latitudes[J]. Hydrographic Surveying and Charting, 2017, 37(4): 19-23
[19]	边少锋, 李厚朴. 大地测量计算机代数分析[M]. 北京: 科学出版社, 2018 Bian Shaofeng, Li Houpu. Computer Algebra Analysis on Geodesy[M]. Beijing: Science Press, 2018
[20]	Bian S F, Chen Y B. Solving an Inverse Problem of a Meridian Arc in Terms of Computer Algebra System[J]. Journal of Surveying Engineering, 2006, 132(1): 7-10 doi: 10.1061/(ASCE)0733-9453(2006)132:1(7)
[21]	丁大正. Mathematica基础与应用[M]. 北京: 电子工业出版社, 2013 Ding Dazheng. Mathematica Foundation and Application[M]. Beijing: Publishing House of Electronics Industry, 2013

施引文献(15)

期刊类型引用(6)

1.	王晗，陈怡霖，季钰姣，杜若琳. 基于多模态唇部状态识别的语音导航抗干扰系统. 江苏大学学报(自然科学版). 2025(01): 82-90 . 百度学术
2.	何进辛，高井祥，潘诚，王永辉. 移动窗口多普勒平滑伪距增强的GNSS/IMU/单目视觉紧组合定位. 测绘通报. 2025(01): 101-106 . 百度学术
3.	王博，武瑞宏，许双安，叶茂，邓川，张卫龙. 低成本亚米级PPP/INS组合导航定位系统研究. 铁道勘察. 2023(02): 8-13+43 . 百度学术
4.	祝浩祈，王甫红，张万威，栾梦杰，程雨欣. 一种城市环境下GNSS/MEMS IMU车载实时精密定位方法. 武汉大学学报(信息科学版). 2023(07): 1232-1240 . 百度学术
5.	朱东琴，王红茹，岳敬轩. 基于EM估计的多模型车载组合导航算法. 仪器仪表学报. 2022(07): 239-246 . 百度学术
6.	李楠，向文豪. 城市环境中无人机作战导航定位研究现状综述. 无人系统技术. 2022(04): 75-87 . 百度学术

其他类型引用(9)

资源附件(0)

表(3)

计量

文章访问数: 586
HTML全文浏览量: 183
PDF下载量: 52
被引次数: 15

1 PPP-RTK/VIO组合模型
1.1 PPP-RTK模型
1.2 INS模型
1.3 VIO模型
1.4 PPP-RTK/VIO半紧组合模型
1.5 模型结构
2 实验与结果分析
2.1 实验说明
2.2 结果分析
3 结语

兰勃脱等角圆锥投影反解不同算法的解析

作者简介: 丁士俊，博士，教授，研究方向为现代测量数据处理理论、测量与地图投影理论与方法。shjding@sgg.whu.edu.cn

通讯作者: 李鹏鹏，博士。276149876@qq.com

计量

出版历程