Low-cost GNSS Heading Determination with Fixed Baseline Constraints
-
摘要: GNSS测向是导航领域提供航向信息的重要手段。当前高性能GNSS测向产品精度可达0.20(°)/m，满足多数应用领域的需求，但其高昂的万元级成本阻碍了大众市场推广。而低成本GNSS接收机的观测量冗余度与质量均弱于高性能GNSS接收机，导致其测向精度与可靠性易受复杂环境影响。为提升低成本GNSS测向性能，本文提出固定基线约束的低成本GNSS测向方法。该方法利用双天线间固定基线约束信息，基于几何基线后验方差-协方差信息构建实时动态检验阈值，可有效增加整周模糊度解算成功率，从而提升低成本GNSS测向的精度与可靠性。通过实测动态数据对所提方法进行测试，与传统方法对比，整周模糊度固定成功率提升了7%，达到了0.40(°)/m的测向精度。Abstract: Objectives: Global Navigation Satellite System (GNSS) heading determination is an important technology in the field of navigation. At present, the accuracy of high-performance GNSS heading determination receivers can reach 0.20(°)/m, which can satisfy the requirements of most applications, but the cost of 10, 000 yuan hinders the mass market promotion. However, the observation redundancy and quality of the low-cost GNSS receivers are weaker than those of the high-performance GNSS receiver, which makes the heading determination accuracy and reliability of the low-cost GNSS receiver vulnerable to complex environments. Methods: In order to improve the performance of low-cost GNSS heading determination, a low-cost GNSS heading determination method with fixed baseline constraints is proposed in this paper. The method deploys fixed baseline constraints between the dual antennas and constructs real-time dynamic detection thresholds based on a posteriori variance-covariance of geometric baseline, which can effectively increase the success rate of integer ambiguity resolution, thus improving the accuracy and reliability of low-cost GNSS heading determination. Results: The proposed method was tested with real kinematic data. Compared to the traditional methods, the success rate of the integer ambiguity resolution was improved by 7%, and the heading determination accuracy of 0.4 (°)/m was achieved. Conclusions: It is shown that the low-cost GNSS vectoring method with fixed baseline constraints proposed in this paper can effectively increase the success rate of whole-period ambiguity solution and improve the accuracy and reliability of low-cost GNSS vectoring.
-
经典平差模型和最小二乘估计理论[1]在大地测量等众多科学研究和工程领域中应用广泛，其中，高斯-马尔科夫模型（Gauss-Markov model, GMM）最为常用，而高斯-赫尔默特模型（GaussHelmert model, GHM）可视为经典平差模型的一般通用形式。在实际应用中，坐标转换、回归模型、数字地面模型和大地测量反演等平差模型的系数矩阵包含随机的观测误差，从而使得GMM扩展为随机系数矩阵的变量含误差（errors-invariables, EIV）模型[2]。文献[3]提出同时顾及观测向量和系数矩阵中随机误差的整体最小二乘（total least squares, TLS）估计算法。TLS的非线性特征导致其受制于计算机技术的发展，直至20世纪80年代，文献[4]将TLS引入数值分析领域并提出奇异值分解算法，TLS才开始广泛应用于各专业领域并取得丰富的研究成果。文献[2]中对TLS进行了改进和扩展；文献[5]从TLS的算法、统计特性和可靠性研究等方面综述了TLS方法的研究进展。当误差相关且精度不等时，采用加权整体最小二乘估计（weighted total least squares, WTLS）方法进行求解。文献[6]研究了基于高斯-牛顿迭代法的WTLS算法，该算法假设权矩阵为特殊情况，得到的解在形式上与最小二乘（least squares, LS）解相同；文献[7]研究了在任意权矩阵的一般情况下的WTLS算法；文献[8]研究了特殊结构下WTLS算法的迭代方法并将其应用于实际场景；文献[9-12]研究了附有等式和不等式约束情况下的WTLS算法；文献[13]研究了稳健WTLS算法。
通过对EIV模型的形式进行变换，文献[14]提出部分EIV(partial EIV, PEIV）模型，提高了系数矩阵仅含部分随机量情况下的计算效率。文献[15]对PEIV模型进行线性化，推导了PEIV模型的LS算法。文献[16-17]从模型的一般性出发，将EIV模型扩展至通用EIV模型，将经典平差的GHM中观测向量的系数矩阵和参数向量的系数矩阵由固定矩阵推广为随机矩阵，涵括随机系数矩阵的各类情况，同时推导了通用EIV模型在任意权矩阵情况下的一般性WTLS算法。
通用EIV模型的非线性使得该算法在估计量较多时计算量大。本文利用非线性平差原理，将通用EIV模型展开后的二阶项纳入平差方程的常数项，从而将其转化为GHM形式，推导出通用EIV模型的线性化整体最小二乘（linearized total least squares, LTLS）算法。相较于WTLS算法，LTLS算法提高了通用EIV模型的计算效率，当参数向量初始值与最优值相差较大时，提升了迭代收敛速度。
1 通用EIV平差模型及其WTLS估计
1.1 通用EIV平差模型
GHM的形式为：
$$ \mathit{\boldsymbol{A}}\left( {\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right) + \mathit{\boldsymbol{BX}} + \mathit{\boldsymbol{w}} = \mathit{\boldsymbol{0}} $$ (1) 式中，y和vy分别为n×1维观测值向量和观测值改正数向量；X为u×1维参数向量；A为观测值向量对应的f×n维系数矩阵；B为参数向量对应的f×u维系数矩阵；w为f×1维常数向量；在经典平差函数模型的定义中，A和B均不含随机误差，为固定矩阵。
当参数向量的系数矩阵B含随机误差时，GHM（式（1））扩展为经典EIV模型。当观测值向量的系数矩阵A和参数向量的系数矩阵B均含随机误差时，GHM（式（1））扩展为通用EIV平差模型[16]:
$$ \left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}} \right)\left( {\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right) + \left( {\mathit{\boldsymbol{B}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}} \right)\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{w}} = \mathit{\boldsymbol{0}} $$ (2) 式中，A和VA分别为观测值向量对应的f×n维系数矩阵及其改正数矩阵；B和VB分别为参数向量对应的f×u维系数矩阵及其改正数矩阵。由于A、B和y均为随机矩阵，则通用EIV的随机模型为：
$$ \mathit{\boldsymbol{L}} = \left[ {\begin{array}{*{20}{c}} {{\mathop{\rm vec}\nolimits} (\mathit{\boldsymbol{A}})}\\ {{\mathop{\rm vec}\nolimits} (\mathit{\boldsymbol{B}})}\\ \mathit{\boldsymbol{y}} \end{array}} \right],\mathit{\boldsymbol{v}} = \left[ {\begin{array}{*{20}{c}} {{\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}} \right)}\\ {{\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}} \right)}\\ {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{A}}}}\\ {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{B}}}}\\ {{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \end{array}} \right] $$ (3) $$ \mathit{\boldsymbol{D}}(\mathit{\boldsymbol{L}}) = \delta _0^2{\mathit{\boldsymbol{P}}^{ - 1}} = \delta _0^2\mathit{\boldsymbol{Q}} = \delta _0^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{A}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{AB}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{Ay}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{BA}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{By}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yA}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yB}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}} \end{array}} \right] $$ (4) 式中，vec(.)表示将矩阵按列向量化；L和v分别为观测数据的k×1维观测值向量及其改正数向量，包括A、B和y中所有观测值及其改正数，其中k=fn+fu+n; P、Q和D(L)分别为L的权矩阵、协因数矩阵和方差协方差矩阵；δ02为单位权方差。
1.2 WTLS估计
根据TLS准则，通用EIV平差模型的求解可转化为最优化估计问题[16]:
$$ \left\{ {\begin{array}{*{20}{l}} {\min {\mathit{\boldsymbol{v}}^{\rm{T}}}\mathit{\boldsymbol{Pv}}}\\ {s.t.{\rm{ }}\left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}} \right)\left( {\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right) + \left( {\mathit{\boldsymbol{B}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}} \right)\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{w}} = \mathit{\boldsymbol{0}}} \end{array}} \right. $$ (5) 相应目标函数为：
$$ \begin{array}{*{20}{c}} {\mathit{\Phi }(\mathit{\boldsymbol{r}},\mathit{\boldsymbol{\lambda }},\mathit{\boldsymbol{X}}) = {\mathit{\boldsymbol{v}}^{\rm{T}}}\mathit{\boldsymbol{Pv}} + 2{\mathit{\boldsymbol{\lambda }}^{\rm{T}}}\left( {\mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} + } \right.}\\ {\left. {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} + \mathit{\boldsymbol{BX}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{X}} + \mathit{\boldsymbol{w}}} \right)} \end{array} $$ (6) 将目标函数对估计量分别求偏导并令其等于0，得到非线性方程组：
$$ \frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{\hat X}}}} = 2\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{\hat \lambda }} + \mathit{\boldsymbol{\hat V}}_\mathit{\boldsymbol{B}}^{\rm{T}}\mathit{\boldsymbol{\hat \lambda }}} \right) = \mathit{\boldsymbol{0}} $$ (7) $$ \frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{\hat v}}}} = 2\left( {\mathit{\boldsymbol{P\hat v}} + {{\mathit{\boldsymbol{\hat C}}}^{\rm{T}}}\mathit{\boldsymbol{\hat \lambda }}} \right) = \mathit{\boldsymbol{0}} $$ (8) $$ \frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{\hat \lambda }}}} = 2(\mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{B\hat X}} + \mathit{\boldsymbol{\hat C\hat v}} + \mathit{\boldsymbol{w}}) = \mathit{\boldsymbol{0}} $$ (9) 式中，$ {\boldsymbol{C}}=\left[{\boldsymbol{y}}^{\mathrm{T}} \otimes {\boldsymbol{I}}_{f} \boldsymbol{X}^{\mathrm{T}} \otimes {\boldsymbol{I}}_{f} A+\boldsymbol{V}_{{\boldsymbol{A}}}\right] ; \hat{{\boldsymbol{v}}}, \hat{\boldsymbol{X}} $分别为观测值向量和参数向量的估计值。
根据式（7）~式（9）可导出：
$$ \mathit{\boldsymbol{\hat v}} = - \mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{\hat C}}^{\rm{T}}}\mathit{\boldsymbol{\hat Q}}_\mathit{\boldsymbol{C}}^{ - 1}(\mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{B\hat X}} + \mathit{\boldsymbol{w}}) $$ (10) $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat X}} = - {{\left[ {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}} + \mathit{\boldsymbol{\hat V}}_\mathit{\boldsymbol{B}}^{\rm{T}}} \right)\mathit{\boldsymbol{\hat Q}}_\mathit{\boldsymbol{C}}^{ - 1}\mathit{\boldsymbol{B}}} \right]}^{ - 1}}\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}} + } \right.}\\ {\left. {\mathit{\boldsymbol{\hat V}}_\mathit{\boldsymbol{B}}^{\rm{T}}} \right)\mathit{\boldsymbol{\hat Q}}_\mathit{\boldsymbol{C}}^{ - 1}(\mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{w}})} \end{array} $$ (11) 式中，$ \hat{\boldsymbol{Q}}_{C}=\hat{\boldsymbol{C}} \boldsymbol{Q} \hat{\boldsymbol{C}}^{\mathrm{T}} $。
以式（2）的LS解作为初始值，根据式（10）和式（11）进行迭代计算可得通用EIV模型的WTLS最优解。
2 通用EIV模型的线性化估计算法
通用EIV模型是非线性模型，观测值矩阵和系数矩阵均为随机量，WTLS算法的计算量随着待估量个数增多将迅速增加。将式（2）展开，利用非线性函数平差原理[18]将二阶项作为模型误差纳入方程的常数项，从而将通用EIV模型转化为线性的GHM，推导出通用EIV模型的LTLS算法。
令X=X0+x，将式（2）展开得：
$$ \begin{array}{c} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{X}}_0} + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right) + \mathit{\boldsymbol{Bx}} + (\mathit{\boldsymbol{w}} + \mathit{\boldsymbol{Ay}} + \\ \left. {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{X}}_0} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{x}}} \right) = \mathit{\boldsymbol{0}} \end{array} $$ (12) 式（12）可表示为：
$$ \begin{array}{*{20}{c}} {\left[ {\left( {{\mathit{\boldsymbol{y}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_\mathit{\boldsymbol{f}}}} \right){\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}} \right) + \left( {\mathit{\boldsymbol{X}}_0^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_\mathit{\boldsymbol{f}}}} \right){\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}} \right) + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}}} \right] + }\\ {\mathit{\boldsymbol{Bx}} + \left( {\mathit{\boldsymbol{w}} + \mathit{\boldsymbol{Ay}} + \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{X}}_0} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{y}}} + {\mathit{\boldsymbol{V}}_\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{x}}} \right) = \mathit{\boldsymbol{0}}} \end{array} $$ (13) 则通用EIV模型的线性化形式为：
$$ {\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{l}}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{l}}} + \mathit{\boldsymbol{Bx}} + {\mathit{\boldsymbol{w}}_\mathit{\boldsymbol{l}}} = \mathit{\boldsymbol{0}} $$ (14) 式中，$ {\boldsymbol{A}}_{{\boldsymbol{l}}}=\left[\begin{array}{llll} {\boldsymbol{y}}^{\mathrm{T}} \otimes {\boldsymbol{I_{f}}} & \boldsymbol{X}_{0}^{\mathrm{T}} & \otimes {\boldsymbol{I_{f}}} & {\boldsymbol{A}} \end{array}\right] ; \ {\boldsymbol{v_{l}}}= \left[\begin{array}{c} \operatorname{vec}\left(\boldsymbol{V}_{{\boldsymbol{A}}}\right) \\ \operatorname{vec}\left(\boldsymbol{V}_{{\boldsymbol{B}}}\right) \\ \boldsymbol{v}_{{\boldsymbol{y}}} \end{array}\right] ; \boldsymbol{w}_{{\boldsymbol{l}}}=\boldsymbol{w}+\boldsymbol{Ay} + {\boldsymbol{B}} \boldsymbol{X}_{0}+\boldsymbol{V}_{{\boldsymbol{A}}} \boldsymbol{v}_{{\boldsymbol{y}}}+\boldsymbol{V}_{{\boldsymbol{B}}} {\boldsymbol{x}}$为常数项，即将式（2）按泰勒展开略去的二阶项VAvy+VBx作为模型误差纳入常数项。
由式（14）可知，线性化后的通用EIV模型与GHM形式一致，可使用最小二乘法得到观测值向量和参数向量：
$$ \mathit{\boldsymbol{\hat x}} = - \mathit{\boldsymbol{\hat N}}_{\mathit{\boldsymbol{bb}}}^{ - 1}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{\hat N}}_{\mathit{\boldsymbol{aa}}}^{ - 1}{\mathit{\boldsymbol{w}}_\mathit{\boldsymbol{l}}} $$ (15) $$ {{\mathit{\boldsymbol{\hat v}}}_\mathit{\boldsymbol{l}}} = - \mathit{\boldsymbol{Q}}{{\mathit{\boldsymbol{\hat A}}}_\mathit{\boldsymbol{l}}}\mathit{\boldsymbol{\hat N}}_{\mathit{\boldsymbol{aa}}}^{ - 1}\left( {\mathit{\boldsymbol{B\hat x}} + {\mathit{\boldsymbol{w}}_\mathit{\boldsymbol{l}}}} \right) $$ (16) 式中，$ \boldsymbol{N}_{{\boldsymbol{a a}}}=\boldsymbol{A}_{{\boldsymbol{l}}} \boldsymbol{Q} \boldsymbol{A}_{{\boldsymbol{l}}}^{\mathrm{T}} ; \boldsymbol{N}_{{\boldsymbol{b b}}}= {\boldsymbol{B}}^{\mathrm{T}} \boldsymbol{N}_{{\boldsymbol{a a}}}^{-1} {\boldsymbol{B}} $; Q为观测值向量的协因数矩阵。
通用EIV模型的LTLS算法步骤如下：
1）将实际模型表示为式（2），将观测数据代入得到A、B和y矩阵，并给出观测值数据的协因数阵Q，包括观测值向量Qy、观测值向量系数矩阵QA、参数向量系数矩阵QB。
2）计算通用EIV模型的LS解作为初始参数解：$ \hat{\boldsymbol{X}}^{0}=-\left({\boldsymbol{B}}^{\mathrm{T}}\left({\boldsymbol{A }}\boldsymbol{Q}_{{\boldsymbol{y}}} {\boldsymbol{A}}^{\mathrm{T}}\right)^{-1} {\boldsymbol{B}}\right)^{-1} {\boldsymbol{B}}^{\mathrm{T}}\left({\boldsymbol{A}} \boldsymbol{Q}_{{\boldsymbol{y}}} {\boldsymbol{A}}^{\mathrm{T}}\right)^{-1} · ({\boldsymbol{A y}}+{\boldsymbol{w}})$, 观测值向量改正数初始值取0。
3）根据式（15）和式（16）进行迭代计算，每次迭代将上一次估计值作为初始值代入新的迭代过程，直至前后两次估计值之差小于设定阈值。
GHM按泰勒级数展开仅包含常数项、一阶项（二阶及以上项全部为零），将二阶项纳入线性化后的常数项，该方法极大地减弱了线性化引起的模型误差。因此，在同样以LS解作为初值的情况下，GHM线性化的LTLS解与WTLS解的收敛性一致。此外，根据文献[19]中EIV模型LS解偏差的研究结果，在当前测量技术手段和观测精度条件下，以有偏的LS解作为初值，能够保证TLS迭代计算收敛，除非出现极特殊情况导致LS初始解严重偏离最优值。
参考GHM，式（2）的LTLS算法估计结果的精度计算式为：
$$ \left\{ {\begin{array}{*{20}{l}} {\hat \delta _0^2 = \mathit{\boldsymbol{v}}_\mathit{\boldsymbol{l}}^{\rm{T}}\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{v}}_\mathit{\boldsymbol{l}}}/\mathit{r}}\\ {\mathit{\boldsymbol{Q}}(\mathit{\boldsymbol{\hat x}}) = {{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{{\left( {{\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{l}}}\mathit{\boldsymbol{QA}}_\mathit{\boldsymbol{l}}^{\rm{T}}} \right)}^{ - 1}}\mathit{\boldsymbol{B}}} \right)}^{ - 1}}}\\ {\mathit{\boldsymbol{Q}}(\mathit{\boldsymbol{\hat L}}) = \mathit{\boldsymbol{Q}} - \mathit{\boldsymbol{QA}}_\mathit{\boldsymbol{l}}^{\rm{T}}\left( {\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{A}}^{ - 1}\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{\hat x}}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{A}}^{ - 1}} \right){\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{l}}}\mathit{\boldsymbol{Q}}} \right)}\\ {\mathit{\boldsymbol{D}}(\mathit{\boldsymbol{\hat x}}) = \hat \delta _0^2\mathit{\boldsymbol{Q}}(\mathit{\boldsymbol{\hat x}})}\\ {\mathit{\boldsymbol{D}}(\mathit{\boldsymbol{\hat L}}) = \hat \delta _0^2\mathit{\boldsymbol{Q}}(\mathit{\boldsymbol{\hat L}})} \end{array}} \right. $$ (17) 式中，多余观测数r=n-t; n为观测值个数；t为必要观测值个数；r与GHM的多余观测数相同；$ \boldsymbol{N}_{{\boldsymbol{A}}}=\boldsymbol{A}_{{\boldsymbol{l}}} \boldsymbol{P}^{-1} \boldsymbol{A}_{{\boldsymbol{l}}}^{\mathrm{T}} $。
3 LTLS算法实例分析
按照LTLS算法步骤设计实验，比较LTLS算法与WTLS算法的计算结果，验证LTLS算法的正确性、高效性和可行性。实验1设计模拟数据，比较分析单组实验结果和1 000组实验统计结果，验证LTLS算法的正确性；实验2在待估计量数目取不同量级时，比较两种算法的计算时间，验证LTLS算法的高效性；实验3通过实例验证LTLS算法的可行性。
3.1 实验1
在通用EIV模型（式（2））中，设置参数真值X=[5 10]T，系数矩阵A和B中随机量的中误差分别设为0.01和0.02，观测向量y的中误差设为0.03。A、B、y和常数向量w的模拟数据如下：
$$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {12.469}&{11.096}&{15.872}&{11.725}\\ {8.883}&{10.291}&{2.929}&{3.666}\\ {12.321}&{1.109}&{6.392}&{15.809}\\ {3.551}&{12.104}&{3.867}&{1.257} \end{array}} \right], $$ $$ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{l}} {10.410}&{17.544}\\ {18.033}&{15.171}\\ {18.631}&{15.855}\\ {15.671}&{11.878} \end{array}} \right], $$ $$ \mathit{\boldsymbol{y}} = \left[ {\begin{array}{*{20}{l}} {27.543}\\ {20.727}\\ {20.839}\\ {25.033} \end{array}} \right],\mathit{\boldsymbol{w}} = \left[ {\begin{array}{*{20}{c}} { - 1425.323}\\ { - 852.619}\\ { - 1142.913}\\ { - 658.407} \end{array}} \right]。 $$ 采用LTLS算法和WTLS算法的估计结果如表 1所示，LTLS参数解与WTLS参数解完全相等，验证了LTLS算法的正确性。
表 1 参数解及其方差估计值Table 1. Parameter Values and Mean Square Deviations算法 参数解 中误差 $ \hat {\boldsymbol{X}} _1 $ $ \hat {\boldsymbol{X}} _2 $ $ σ _{\hat{\boldsymbol{X}}_1} $ $ σ_ {\hat{\boldsymbol{X}}_2} $ LTLS算法 5.012 551 9.994 964 0.039 9 0.051 9 WTLS算法 5.012 551 9.994 964 0.039 9 0.051 9 为了进一步验证LTLS算法的正确性，首先采用模拟的1 000组数据计算LTLS和WTLS参数解的均值$ {\rm{avg}} (\hat {\boldsymbol{X}} _1) 和 {\rm{avg}} (\hat {\boldsymbol{X}} _2) $，并将参数解均值代入式（17），求得LTLS参数解的协因数阵估值$ {\boldsymbol{Q}}(\hat {\boldsymbol{X}} )$，然后利用参数真值求得精确的参数协因数阵$ \overline{\boldsymbol{Q}}(\hat{\boldsymbol{X}})={\boldsymbol{E}}_{{\boldsymbol{X}}}^{\mathrm{T}} \boldsymbol{E}_{{\boldsymbol{X}}} /(m-n) , {\boldsymbol{E}}_{{\boldsymbol{X}}}=\left[\begin{array}{ll} \hat{{\boldsymbol{X}}}_{1}^{(j)}-5 & \hat{{\boldsymbol{X}}}_{2}^{(j)}-10 \end{array}\right] $，计算结果见表 2。1 000组数据的LTLS参数解和WTLS参数解与参数真值偏差的统计分析见图 1和图 2。
表 2 1000组实验的参数解均值和协因数阵Table 2. Average Parameter Values and Co⁃variance Matrix in 1 000 Experiments算法 参数解均值 协因数阵 $ {\rm{avg}} ( \hat{\boldsymbol{X}}_1) $ $ {\rm{avg}} ( \hat{\boldsymbol{X}}_2) $ $ {\boldsymbol{Q }}( \hat{\boldsymbol{X}} ) $ $ \bar{ \boldsymbol{Q }}( \hat{\boldsymbol{X}} ) $ LTLS算法 5.009 427 9.998 257 0.004 3 −0.005 1 0.004 2 −0.004 9 −0.005 1 0.007 3 −0.004 9 0.007 2 WTLS算法 5.009 427 9.998 257 0.004 3 −0.005 1 0.004 2 −0.004 9 −0.005 1 0.007 3 −0.004 9 0.007 2 ![]()
图 1 1 000组实验数据的LTLS参数解与参数真值偏差统计图Figure 1. Statistical Graph of Deviations Between Truth Values and Parameter Values Soluted by LTLS in 1 000 Experiments![]()
图 2 1 000组实验数据的WTLS参数解与参数真值偏差统计图Figure 2. Statistical Graph of Deviations Between Truth Values and Parameter Values Soluted by WTLS in 1 000 Experiments由表 2、图 1和图 2可以看出，LTLS算法与WTLS算法的统计结果完全一致，说明在每次实验中两种算法所求参数解均一致，验证了LTLS算法的正确性。同时两种计算协因数阵的方法结果非常相近，验证了参数协因数阵一阶近似估计公式（17）的有效性[17]。
3.2 实验2
为分析LTLS算法的计算效率，设计通用EIV模型中待估计量个数在不同的数量级情况，采用LTLS算法和WTLS算法计算100组模拟数据的平均迭代次数N、平均解算时间t和减少比例（LTLS算法较WTLS算法减少的平均解算时间与WTLS算法平均解算时间之比），结果见表 3。
表 3 LTLS算法和WTLS算法计算效率的比较Table 3. Comparison of Computational Efficiency Between LTLS Algorithm and WTLS Algorithm待估量数量 NLTLS NWTLS tLTLS tWTLS 减少比例/% 10 5.18 4.5 0.224 ms 0.191 ms − 100 5.04 6.36 0.681 ms 0.724 ms 5.9 1 000 5 6.75 0.033 s 0.044 s 25.0 10 000 5 7.5 2.681 s 3.909 s 31.4 从表 3可以看出，两种算法每次迭代的平均时间基本一致；当模型估计量数量较少时，两种算法的效率基本相当，随着估计量的数量级逐渐增大，LTLS算法的效率高于WTLS算法。原因在于GHM按泰勒级数展开后，仅包含常数项、一阶项（二阶及以上项全部为零），LTLS算法将二阶项纳入常数项，减小了线性化引起的模型误差，迭代计算收敛更快，迭代次数减少，使得计算效率提高。
3.3 实验3
本文采用的摄影测量实例示意图如图 3所示，由3个地面摄像机S1、S2和S3拍摄两个目标点P1和P2组成，相机主距f=100 mm，距离l1、l2、l3、l4、l5、l6、y1、y2的观测值和中误差见表 4。
表 4 距离观测值及其中误差Table 4. Distance Observations and Standard Deviations统计项 l1/mm l2/mm l3/mm l4/mm l5/mm l6/mm y1/m y2/m 观测值 14.1 16.6 6.1 7.1 22.1 26.3 10.0 8.0 中误差 0.10 0.10 0.10 0.10 0.10 0.10 0.05 0.05 根据图 3可得到误差方程：
$$ \left\{\begin{array}{l} l_{1} x_{2}-f x_{1}=0 \\ l_{2} x_{4}-f x_{3}=0 \\ l_{3} x_{2}+f y_{1}+f x_{1}=0 \\ l_{4} x_{4}+f y_{1}+f x_{3}=0 \\ l_{5} x_{2}-f y_{1}-f y_{2}+f x_{1}=0 \\ l_{6} x_{4}-f y_{1}-f y_{2}+f x_{3}=0 \end{array}\right. $$ (18) 构建通用EIV模型来估计点P1和点P2的坐标，由误差方程可得：
$$ \begin{array}{l} \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} 0&0\\ 0&0\\ { - f}&0\\ { - f}&0\\ { - f}&{ - f}\\ { - f}&{ - f} \end{array}} \right],\mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} { - f}&{{l_1}}&0&0\\ 0&0&{ - f}&{{l_2}}\\ f&{{l_3}}&0&0\\ 0&0&f&{{l_4}}\\ f&{{l_5}}&0&0\\ 0&0&f&{{l_6}} \end{array}} \right],\\ \mathit{\boldsymbol{y}} = \left[ {\begin{array}{*{20}{l}} {{y_1}}\\ {{y_2}} \end{array}} \right],\mathit{\boldsymbol{x}} = \left[ {\begin{array}{*{20}{l}} {{x_1}}\\ {{x_2}}\\ {{x_3}}\\ {{x_4}} \end{array}} \right],\mathit{\boldsymbol{w}} = \left[ {\begin{array}{*{20}{l}} 0\\ 0\\ 0\\ 0\\ 0\\ 0 \end{array}} \right] \end{array} $$ (19) 由式（19）可知矩阵A不含有随机误差，则该矩阵的改正数矩阵为零矩阵。采用LTLS算法和WTLS算法解算，设置阈值为1×10-8，结果见表 5和表 6。两种方法所得坐标估值和距离观测值估值完全一致，该实例表示为通用EIV模型时系数矩阵中待估计量较少，所以两个算法收敛速度相差不大。
表 5 点P1和点P2的坐标估值/mTable 5. Coordinate Estimates of P1 and P2 /m算法 坐标估值 $ \hat x_1 $ $ \hat x_2 $ $ \hat x_3 $ $ \hat x_4 $ LTLS算法 6.995 056 5 49.715 632 6.981 465 5 41.968 315 9 WTLS算法 6.995 056 5 49.715 632 6.981 465 5 41.968 315 9 表 6 距离观测值估值Table 6. Estimation of Distance Observations算法 $ \hat l_1 $/mm $ \hat l_2 $/mm $ \hat l_3 $/mm $ \hat l_4 $/mm $ \hat l_5 $5/mm $ \hat l_6 $/mm $ \hat y_1 $/m $ \hat y_2 $/m LTLS算法 14.070 1 16.635 1 6.032 4 7.178 4 22.137 7 26.256 7 9.994 1 8.006 8 WTLS算法 14.070 1 16.635 1 6.032 4 7.178 4 22.137 7 26.256 7 9.994 1 8.006 8 4 结语
本文将通用EIV函数模型展开后的二阶项纳入模型的常数项，将通用EIV模型表示为线性形式的GHM，推导出通用EIV模型的线性化整体最小二乘算法和近似精度估计公式。实验结果表明，通用EIV模型的LTLS算法与WTLS算法结果一致，验证了该算法的正确性。此外，LTLS算法估计精度公式和WTLS估计精度公式均为一阶近似精度，因此两种算法参数的估计精度相同。当通用EIV模型的待估量数量较多时，LTLS算法比WTLS计算效率更高，在处理海量数据时更具有优势。
-
[1] Sun R, Cheng Q, Wang J. Precise Vehicle Dynamic Heading and Pitch Angle Estimation Using Time-differenced Measurements from a Single GNSS Antenna[J]. GPS Solutions, 2020, 24(3):1-9
[2] Cai Xiaobo, Xu Houze, Wang Yong, et al. Direct Attitude Determination Method Based on Vehicle-mounted Three-antenna GNSS and the Accuracy Evaluation[J]. Geomatics and Information Science of Wuhan University, 2018, 43(6):820-825(蔡小波, 许厚泽, 王勇, 等. 车载三天线GNSS的直接法定姿及精度评估[J]. 武汉大学学报(信息科学版), 2018, 43(6):820-825) [3] Ding W, Sun W, Gao Y, et al. Carrier Phase-Based Precise Heading and Pitch Estimation Using a Low-Cost GNSS Receiver[J]. Remote Sensing, 2021, 13(18):3642
[4] Teunissen P, Giorgi G, Buist P. Testing of a New Single-frequency GNSS Carrier Phase Attitude Determination Method:Land, Ship, and Aircraft Experiments[J]. GPS Solutions, 2011, 15(1):15-28
[5] Li N, Zhao L, Li L, et al. Integrity Monitoring of High-accuracy GNSS-based Attitude Determination[J]. GPS Solutions, 2018, 22(4):1-3
[6] Wang X, Yao Y, Xu C, et al. An Improved Single-epoch Attitude Determination Method for Low-cost Single-frequency GNSS Receivers[J]. Remote Sensing, 2021, 13(14):2746
[7] Odolinski R, Teunissen P. Single-frequency, Dual-GNSS Versus Dual-frequency, Single-GNSS:a Low-cost and High-grade Receivers GPS-BDS RTK Analysis[J]. Journal of Geodesy, 2016, 90(11):1255-1278
[8] Lachapelle G, Gratton P. GNSS Precise Point Positioning with Android Smartphones and Comparison with High Performance Receivers[C]//Proc. of the 2019 IEEE International Conference on Signal, Information and Data Processing (ICSIDP). Chongqing, China:IEEE, 2019:1-9.
[9] Zhang Xiaohong, Tao Xianlu, Wang Yingzhe, et al. MEMS-Enhanced Smartphone GNSS High-Precision Positioning for Vehicular Navigation in Urban Conditions[J]. Geomatics and Information Science of Wuhan University, 2022, 47(10):1740-1749(张小红, 陶贤露, 王颖喆, 等. 城市场景智能手机GNSS/MEMS融合车载高精度定位[J]. 武汉大学学报(信息科学版), 2022, 47(10):1740-1749) [10] Yang Y, Mao X, Tian W. A Novel Method for Low-cost MIMU Aiding GNSS Attitude Determination[J]. Measurement Science and Technology, 2016, 27(7):075003
[11] Rychlicki M, Kasprzyk Z, Rosiński A. Analysis of Accuracy and Reliability of Different Types of GPS Receivers[J]. Sensors, 2020, 20(22):1-14
[12] Zhu F, Hu Z, Liu W, et al. Dual-antenna GNSS Integrated with MEMS for Reliable and Continuous Attitude Determination in Challenged Environments[J]. IEEE Sensors Journal, 2019, 19(9):3449-3461
[13] Tang Weiming, Sun Hongxing, Liu Jingnan. Ambiguity Resolution of Single Epoch Single Frequency Data with Baseline Length Constraint Using LAMBDA Algorithm[J]. Geomatics and Information Science of Wuhan University, 2005, 30(5):444-446(唐卫明, 孙红星, 刘经南. 附有基线长度约束的单频数据单历元LAMBDA方法整周模糊度确定[J]. 武汉大学学报(信息科学版), 2005, 30(5):444-446) [14] Ma L, Lu L, Zhu F, et al. Baseline length Constraint Approaches for Enhancing GNSS Ambiguity Resolution:Comparative Study[J]. GPS Solutions, 2021, 25(2):1-5
[15] Wu S, Zhao X, Zhang L, et al. Improving Reliability and Efficiency of RTK Ambiguity Resolution with Reference Antenna Array:BDS+ GPS Analysis and Test[J]. Journal of Geodesy. 2019, 93(9):1297-1311
[16] Tang Jiaming, Chai Yanju, Wen Debao, et al. Algorithm of Simultaneous Ambiguity Resolution of GNSS Multi-Antenna Network for Single Epoch[J]. Journal of Geodesy and Geodynamics, 2019, 39(03):262-268(汤佳明, 柴艳菊, 闻德保, 等. GNSS多天线基线网单历元模糊度同步解算法[J]. 大地测量与地球动力学, 2019, 39(03):262-268) [17] Nie Zhixi, Wang Zhenjie, Ou Jikun, et al. On the Effect of Linearization and Approximation of Nonlinear Baseline Length Constraint for Ambiguity Resolution[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(2):168-173(聂志喜, 王振杰, 欧吉坤, 等. 非线性基线长约束条件线性化近似对模糊度解算的影响[J]. 测绘学报, 2015, 44(02):168-173) [18] Li Y, Zhao L, Jia C, et al. Low-Cost Dual-Antenna GNSS Precision Heading Determination Method with Baseline Length Constraint[C]//Proc. of the China Satellite Navigation Conference (CSNC 2021) Proceedings. Nanchang, China:Springer, 2021:389-402
[19] Teunissen P. The LAMBDA Method for the GNSS Compass[J]. Artificial Satellites, 2006, 41(3):89-103
[20] Zhou Xiaoqing, Shan Hongyi, Lu Liguo, et al. Performance Comparison Between Two Baseline Length Constraint Methods of Ambiguity Resolution[J]. Science of Surveying and Mapping, 2016, 41(6):53-58(周晓青, 单弘煜, 卢立果, 等. 附有长度约束的模糊度解算方法比较[J]. 测绘科学, 2016, 41(6):53-8.) [21] Jia Chun, Zhao Lin, Li Liang, et al. BDS Triple-frequency Tightly Coupled Short-baseline RTK Method by Calibrating the Between-receiver Inter-frequency Biases[J]. Scientia Sinica Terrae, 2020, 50(1):90-103(贾春, 赵琳, 李亮, 等. 改正接收机频间偏差的短基线北斗三频紧组合RTK方法[J]. 中国科学:地球科学. 2020, 50(1):90-103) [22] Liu Jingnan, Deng Chenlong, Tang Weiming. Review of GNSS Ambiguity Validation Theory[J]. Geomatics and Information Science of Wuhan University, 2014, 39(9):1009-1016(刘经南, 邓辰龙, 唐卫明. GNSS整周模糊度确认理论方法研究进展[J]. 武汉大学学报(信息科学版), 2014, 39(9):1009-1016) -
期刊类型引用(6)
1. 杨根新，王友昆，谢正明. 基于广义EIV模型的矿区高程异常的无缝推估算法. 工程勘察. 2023(08): 46-51 . 
百度学术
2. Jianjun ZHU，Leyang WANG，Jun HU，Bofeng LI，Haiqiang FU，Yibin YAO. Recent Advances in the Geodesy Data Processing. Journal of Geodesy and Geoinformation Science. 2023(03): 33-45 . 必应学术
3. 翁烨，陈丽，王岩. 线性化通用EIV平差模型的正则化解法. 勘察科学技术. 2023(05): 1-5 . 百度学术
4. 戴中东，孟良，高永攀，项伟. 加权整体最小二乘坐标匹配算法在机场道面测量中的应用. 测绘地理信息. 2022(02): 61-66 . 百度学术
5. 翁烨，邵德盛，甘淑. 线性化通用EIV平差模型的岭估计解法. 全球定位系统. 2022(02): 82-89 . 百度学术
6. 翁烨，邵德盛. 病态加权总体最小二乘的广义岭估计解法. 全球定位系统. 2021(06): 84-89 . 百度学术
其他类型引用(2)
下载:
计量
- 文章访问数: 394
- HTML全文浏览量: 22
- PDF下载量: 41
- 被引次数: 8
