留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

测距定位方程参数估计的Frozen-Barycentre算法

孙振 曲国庆 苏晓庆 杜存鹏 邓晓景

孙振, 曲国庆, 苏晓庆, 杜存鹏, 邓晓景. 测距定位方程参数估计的Frozen-Barycentre算法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(9): 1478-1484. doi: 10.13203/j.whugis20180129
引用本文: 孙振, 曲国庆, 苏晓庆, 杜存鹏, 邓晓景. 测距定位方程参数估计的Frozen-Barycentre算法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(9): 1478-1484. doi: 10.13203/j.whugis20180129
SUN Zhen, QU Guoqing, SU Xiaoqing, DU Cunpeng, DENG Xiaojing. Frozen-Barycentre Algorithm for Solving Distance Equations[J]. Geomatics and Information Science of Wuhan University, 2020, 45(9): 1478-1484. doi: 10.13203/j.whugis20180129
Citation: SUN Zhen, QU Guoqing, SU Xiaoqing, DU Cunpeng, DENG Xiaojing. Frozen-Barycentre Algorithm for Solving Distance Equations[J]. Geomatics and Information Science of Wuhan University, 2020, 45(9): 1478-1484. doi: 10.13203/j.whugis20180129

测距定位方程参数估计的Frozen-Barycentre算法

doi: 10.13203/j.whugis20180129
基金项目: 

国家高技术研究发展计划 2016YDB051700

国家自然科学基金 41674014

国家自然科学基金 41704003

山东省自然科学基金 BS2014HZ016

详细信息
    作者简介:

    孙振, 硕士生, 主要从事现代测量数据处理研究。98625410@qq.com

    通讯作者: 曲国庆, 博士, 教授。qgq@sdut.edu.cn
  • 中图分类号: P207

Frozen-Barycentre Algorithm for Solving Distance Equations

Funds: 

The National High-tech R & D Program of China 2016YDB051700

the National Natural Science Foundation of China 41674014

the National Natural Science Foundation of China 41704003

the Natural Science Foundation of Shandong Province BS2014HZ016

More Information
    Author Bio:

    SUN Zhen, postgraduate, specializes in the modern measurement data processing. E-mail:98625410@qq.com

    Corresponding author: QU Guoqing, PhD, professor. E-mail: qgq@sdut.edu.cn
  • 摘要: 针对传统病态非线性最小二乘求解不稳定且可靠性低的特点, 基于测距定位方程最小二乘解性质, 提出了一种Frozen-Barycentre迭代法。该方法将萨玛斯基应用于重心迭代法, 实现了内迭代和外迭代的转换, 通过减少导数计算量节省运算时间, 提高重心迭代法的收敛效率。并采用模拟数据和水下定位实测数据, 验证了该方法的数值收敛解优于线性化平差估计解, 收敛效率优于重心迭代法。
  • 图  1  高斯—牛顿法的点位迭代序列

    Figure  1.  The Sequence of Gauss-Newton Method

    图  2  Barycentre和Frozen-Barycentre法的点位迭代序列

    Figure  2.  The Sequences of Barycentre and Frozen-Barycentre Methods

    图  3  测量船围绕应答器实测航线

    Figure  3.  The Ship's Measured Route Around the Transponder

    图  4  15个相邻数据的已知点坐标

    Figure  4.  Known Point Coordinates for 15 Adjacent Data

    图  5  Barycentre和Frozen-Barycentre法的点位迭代序列

    Figure  5.  The Sequences of Barycentre and Frozen-Barycentre Methods

    表  1  短程距离观测控制点坐标和观测距离/m

    Table  1.   Trilateration Network Known Coordinates and Observation Distances/m

    点号 坐标 距离
    X Y Z
    P1 10.000 20.000 6.000 22.180
    P2 8.000 30.000 0.000 30.960
    P3 10.000 -20.000 6.000 23.170
    P4 8.000 -30.000 0.000 30.880
    P5 -10.000 -20.000 -6.000 23.070
    P6 -8.000 -30.000 0.000 30.440
    P7 -10.000 20.000 -6.000 22.700
    P8 -8.000 30.000 0.000 30.720
    下载: 导出CSV

    表  2  不同方法的解算结果

    Table  2.   Calculation Results of Different Methods

    方法 t 数值收敛解$\mathit{\boldsymbol{\hat x}}$/m 时间/s 外循环k 残差$\left\| {\Delta \mathit{\boldsymbol{\hat x}}} \right\|/{\bf{m}}$
    线性化平差 -0.648 7 0.1237 1.4791 1.6198
    Barycentre -0.214 1 0.1238 0.5493 0.130909 987 0.6024
    Frozen-Barycentre 2 -0.214 1 0.1238 0.5493 0.119818 491 0.6024
    3 -0.214 1 0.1238 0.5493 0.097318 326 0.6024
    4 -0.214 1 0.1238 0.5493 0.085686 243 0.6024
    5 -0.214 1 0.1238 0.5493 0.078891 193 0.6024
    6 -0.214 1 0.1238 0.5493 0.073928 160 0.6024
    7 -0.214 1 0.1238 0.5493 0.073005 137 0.6024
    8 -0.214 1 0.1238 0.5493 0.067277 119 0.6024
    9 -0.214 1 0.1238 0.5493 0.065978 105 0.6024
    10 -0.214 1 0.1238 0.5493 0.064174 94 0.6024
    下载: 导出CSV

    表  3  不同算法的水下定位解算结果

    Table  3.   Results of Convergence of Barycentre and Frozen-Barycentre Algorithms

    方法 t 数值收敛解/m 时间/s k 残差$\left\| {\Delta \mathit{\boldsymbol{\hat x}}} \right\|/{\rm{m}}$
    线性化平差 2 553993.560 3 -145 334.776 5 371795 350.815 7 371793 904.189
    Barycentre 2 438949.021 9 492011.361 4 2142.926 3 74.9428 25 032 135.102
    Frozen-Barycentre 2 2 438949.021 9 492011.361 4 2142.926 3 20.2227 12 518 135.102
    3 2 438949.021 9 492011.361 4 2142.926 3 10.3040 8 347 135.102
    4 2 438949.021 9 492011.361 4 2142.926 3 7.2169 6 262 135.102
    5 2 438949.021 9 492011.361 4 2142.926 3 5.0740 5 010 135.102
    6 2 438949.021 9 492011.361 4 2142.926 3 4.3751 4 176 135.102
    7 2 438949.021 9 492011.361 4 2142.926 3 3.5661 3 580 135.102
    8 2 438949.021 9 492011.361 4 2142.926 3 3.2884 3 133 135.102
    9 2 438949.021 9 492011.361 4 2142.926 3 2.9414 2 786 135.102
    10 2 438949.021 9 492011.361 4 2142.926 3 2.7825 2 508 135.102
    下载: 导出CSV
  • [1] 王新洲.非线性模型线性近似的容许曲率[J].武汉大学学报·信息科学版, 1997, 22(2):119-121 http://ch.whu.edu.cn/article/id/4176

    Wang Xinzhou.Linear Approximation of Non-linear Models Allow Curvature[J].Geomatics and Information Science of Wuhan University, 1997, 22(2):119-121 http://ch.whu.edu.cn/article/id/4176
    [2] 王志忠, 阳可奇.非线性模型的线性近似条件研究[J].中南大学学报(自然科学版), 1999(3):230-233 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK199901143619

    Wang Zhizhong, Yang Keqi. Study on Non-linear Model of Linear Approximation[J]. Journal of Central South University:Science and Technology, 1999(3):230-233 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK199901143619
    [3] Mao X, Wada M, Hashimoto H. Nonlinear Iterative Algorithm for GPS Positioning with Bias Model[C].The International IEEE Conference on Intelligent Transportation Systems, Washington, USA, 2004
    [4] Sirola N. Closed-Form Algorithms in Mobile Positioning: Myths and Misconceptions[C]. 7th Workshop on IEEE Positioning Navigation and Communication, Dresden, Germany, 2010
    [5] Gerritsma D I M I. Development and Analysis of Least-Squares Spectral Element Methods for Non-linear Hyperbolic Problems[J]. American Journal of Psychology, 1995, 108, DOI: 10.2307/1423127
    [6] Grafarend E, Schaffrin B. The Geometry of Nonlinear Adjustment—The Planar Trisection Problem[J].Department of Geodesy, 1989:149-172
    [7] Takahama T, Sakai S, Iwane N. Constrained Optimization by the ε, Constrained Hybrid Algorithm of Particle Swarm Optimization and Genetic Algorithm[M]// AI 2005: Advances in Artificial Intelligence. Berlin, Heidelberg :Springer, 2005:389-400
    [8] Teunissen P. Applications of Linear and Nonlinear Models: Fixed Effects, Random Effects, and Total Least Squares[J]. Surveyor, 2012, 58(2):339-340
    [9] Kaltenbacher B, Hofmann B. Convergence Rates for the Iteratively Regularized Gauss-Newton Method in Banach Spaces[J]. Inverse Problems, 2010, 26(3):35 007-35 021 doi:  10.1088/0266-5611/26/3/035007
    [10] Noftsger S, St-Pierre N R. Closed-form Algorithms in Mobile Positioning: Myths and Misconceptions[C]. IEEE Positioning Navigation and Communication, Dresden, Germany, 2010
    [11] 薛树强, 杨元喜, 党亚民.测距定位方程非线性平差的封闭牛顿迭代公式[J].测绘学报, 2014, 43(8):771-777

    Xue Shuqiang, Yang Yuanxi, Dang Yamin. A Closed-form of Newton Iterative Formula for Nonlinear Adjustment of Distance Equations[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(8):771-777
    [12] Xue S, Dang Y, Liu J, et al. Bias Estimation and Correction for Triangle-Based Surface Area Calculations[J]. International Journal of Geographical Information Science, 2016, 30(11):2 155-2 170 doi:  10.1080/13658816.2016.1162795
    [13] Xue Shuqiang, Yang Yuanxi, Dang Yamin. Barycentre Method for Solving Distance Equations[J]. Survey Review, 2016, 48(348):188-194 doi:  10.1179/1752270615Y.0000000020
    [14] Xue S Q, Yang Y X. Gaussâ-Jacobi Combinatorial Adjustment and Its Modification[J]. Empire Survey Review, 2014, 46(337):298-304 doi:  10.1179/1752270613Y.0000000086
    [15] Amat S, Argyros I, Busquier S, et al. On Two High Order Families of Frozen Newton Type Methods[J]. Numerical Linear Algebra with Applications, 2017(107):e2126 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=10.1002/nla.2126
    [16] Xu J, Han B, Li L. Frozen Landweber Iteration for Nonlinear Ill-Posed Problems[J]. Acta Mathematicae Applicatae Sinica, 2007, 23(2):329-336 doi:  10.1007/s10255-007-0375-2
  • 加载中
图(5) / 表(3)
计量
  • 文章访问数:  39
  • HTML全文浏览量:  11
  • PDF下载量:  23
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-08-18
  • 刊出日期:  2020-09-05

测距定位方程参数估计的Frozen-Barycentre算法

doi: 10.13203/j.whugis20180129
    基金项目:

    国家高技术研究发展计划 2016YDB051700

    国家自然科学基金 41674014

    国家自然科学基金 41704003

    山东省自然科学基金 BS2014HZ016

    作者简介:

    孙振, 硕士生, 主要从事现代测量数据处理研究。98625410@qq.com

    通讯作者: 曲国庆, 博士, 教授。qgq@sdut.edu.cn
  • 中图分类号: P207

摘要: 针对传统病态非线性最小二乘求解不稳定且可靠性低的特点, 基于测距定位方程最小二乘解性质, 提出了一种Frozen-Barycentre迭代法。该方法将萨玛斯基应用于重心迭代法, 实现了内迭代和外迭代的转换, 通过减少导数计算量节省运算时间, 提高重心迭代法的收敛效率。并采用模拟数据和水下定位实测数据, 验证了该方法的数值收敛解优于线性化平差估计解, 收敛效率优于重心迭代法。

English Abstract

孙振, 曲国庆, 苏晓庆, 杜存鹏, 邓晓景. 测距定位方程参数估计的Frozen-Barycentre算法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(9): 1478-1484. doi: 10.13203/j.whugis20180129
引用本文: 孙振, 曲国庆, 苏晓庆, 杜存鹏, 邓晓景. 测距定位方程参数估计的Frozen-Barycentre算法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(9): 1478-1484. doi: 10.13203/j.whugis20180129
SUN Zhen, QU Guoqing, SU Xiaoqing, DU Cunpeng, DENG Xiaojing. Frozen-Barycentre Algorithm for Solving Distance Equations[J]. Geomatics and Information Science of Wuhan University, 2020, 45(9): 1478-1484. doi: 10.13203/j.whugis20180129
Citation: SUN Zhen, QU Guoqing, SU Xiaoqing, DU Cunpeng, DENG Xiaojing. Frozen-Barycentre Algorithm for Solving Distance Equations[J]. Geomatics and Information Science of Wuhan University, 2020, 45(9): 1478-1484. doi: 10.13203/j.whugis20180129
  • 由位置信息构成的距离观测方程为非线性方程, 多采用迭代方法进行计算, 这就要求数值算法能够快速、稳定收敛到预期解。研究表明, 非线性扰动主要来源于线性近似时系数矩阵的扰动、附加的截断误差及正交过程, 初始值的选取和近似正交都可能导致线性近似引起的系数矩阵扰动和截断误差增大超出容许范围[1-2]。特别是在非线性模型病态时, 迭代计算不稳定会导致平差结果失真。根据问题的病态和秩亏程度选择合适的平差模型显得尤为重要。为解决非线性最小二乘不适定问题, 有关学者提出了一系列数值迭代方法, 如阻尼最小二乘法、非线性最小二乘同伦法和正则化牛顿法等, 目的是加速非线性最小二乘迭代速度、避免迭代矩阵求逆或者改善迭代计算稳定性[3-5]

    大地测量与卫星导航定位中, 需要求解距离观测方程以获得位置参数, 该方程为超定非线性方程, 需采用非线性最小二乘方法进行参数估计[6-7]。虽然牛顿法能够处理非线性最小二乘问题, 但由于其储存成本等因素, 高斯牛顿迭代法[8-9]一般更为常用。若距离观测方程存在不适定问题时, 高斯牛顿法会由于观测数据误差扰动而出现强烈的不稳定特征, 甚至无法收敛而解算失效, 该问题在移动测量定位领域中已进行探讨[10]。文献[11]联合牛顿法和高斯—牛顿法, 提出了一种封闭式牛顿法来处理距离观测方程。可是, 在处理病态问题时, 牛顿型数值算法过多依赖于初始点精度, 会由于矩阵求逆的奇异性和计算舍入误差而无法收敛, 以致解算失效[12]

    对于测距定位问题, 非线性最小二乘解是观测向量末端以观测权为质量质点系的重心, 基于该条件发展了简单且无需矩阵求逆的稳定数值方法, 即重心迭代法(Barycentre)[13], 然而, 该方法收敛速度可能非常缓慢, 具有较低的收敛效率。为此, 本文提出了一种Frozen-Barycentre迭代法来处理距离观测方程, 提高Barycentre迭代法的收敛效率。

    • 距离观测模型在导航定位中具有重要的地位, 其观测方程为:

      $${L_i} = {d_i}\left( \mathit{\boldsymbol{x}} \right) + {\varepsilon _i}$$ (1)

      式中, ${d_i}\left( \mathit{\boldsymbol{x}} \right) = \left\| {{x_i} - {x_2}} \right\| = \sqrt {\mathop \sum \limits_j^m {{({x_j} - {x_{ij}})}^2}} $为欧氏距离; $\mathit{\boldsymbol{x}} = \left[ {{x_1}{\rm{}}{x_2}{\rm{}} \cdots {\rm{}}{x_m}} \right] \in {R^m}$为待测点构成的向量; ${\varepsilon _i}$为观测误差。

      非线性最小二乘平差实质是求泛函极小值的最优化问题, 常用来处理式(1), 则非线性最小二乘参数估计解为:

      $$\mathit{\boldsymbol{\hat x}} = {\rm{argmin}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( \mathit{\boldsymbol{x}} \right)$$ (2)

      式中,

      $$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( \mathit{\boldsymbol{x}} \right) = \frac{1}{2}{\mathit{\boldsymbol{V}}^{\rm{T}}}\left( \mathit{\boldsymbol{x}} \right)\mathit{\boldsymbol{PV}}\left( \mathit{\boldsymbol{x}} \right)$$ (3)

      其中, $\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}:{R^n} \mapsto R$为非线性最小二乘目标函数; $\mathit{\boldsymbol{P}} = {\rm{diag}}\left( {{p_1}{\rm{}}{p_2}{\rm{}} \cdots {\rm{}}{p_n}} \right)$为观测权阵; 残差向量V为:

      $$\mathit{\boldsymbol{V}}\left( \mathit{\boldsymbol{x}} \right) = \mathit{\boldsymbol{d}}\left( \mathit{\boldsymbol{x}} \right) - \mathit{\boldsymbol{L}}$$ (4)

      运用无约束非线性最小二乘求解式(1), 非线性最小二乘目标函数(式(3))满足如下正交条件:

      $$\mathit{\boldsymbol{h}}\left( {\mathit{\boldsymbol{\hat x}}} \right) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} '}}\left( {\mathit{\boldsymbol{\hat x}}} \right) = {\mathit{\boldsymbol{J}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\hat x}}} \right)\mathit{\boldsymbol{PV}}\left( {\mathit{\boldsymbol{\hat x}}} \right) = 0$$ (5)

      式中, $\mathit{\boldsymbol{J}}\left( {\mathit{\boldsymbol{\hat x}}} \right)$为雅克比矩阵:

      $$\mathit{\boldsymbol{J}}\left( {\mathit{\boldsymbol{\hat x}}} \right) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_1}}\\ {{\mathit{\boldsymbol{e}}_2}}\\ \vdots \\ {{\mathit{\boldsymbol{e}}_n}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{e_{11}}}&{{e_{12}}}& \cdots &{{e_{1m}}}\\ {{e_{21}}}&{{e_{22}}}& \cdots &{{e_{2m}}}\\ \vdots & \vdots &{}& \vdots \\ {{e_{n1}}}&{{e_{n2}}}& \cdots &{{e_{nm}}} \end{array}} \right]$$ (6)

      其中, ${e_{ij}} = \frac{{{x_{ij}} - {{\hat x}_j}}}{{\hat x - {x_i}_2}}$为矢量${\mathit{\boldsymbol{e}}_i} = \frac{{{\mathit{\boldsymbol{x}}_i} - \mathit{\boldsymbol{\hat x}}}}{{d\left( {\mathit{\boldsymbol{\hat x}}} \right)}}$第j个观测方向余弦。由于$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( \mathit{\boldsymbol{x}} \right)$连续可微, 可采用最速下降法来处理正交条件方程(5)。最速下降法的迭代形式为[14]

      $${\mathit{\boldsymbol{\hat x}}_{k + 1}} = {\mathit{\boldsymbol{\hat x}}_k} + \gamma \left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)\mathit{\boldsymbol{G}}\left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)\mathit{\boldsymbol{\psi }}\left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)$$ (7)

      式中, $\gamma \left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)$为线性搜索函数; $\mathit{\boldsymbol{G}}\left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)$为正定矩阵; $\mathit{\boldsymbol{\psi }}\left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right) = \mathit{\boldsymbol{h}}\left( {{\mathit{\boldsymbol{x}}_k}} \right)$。应用如下不动点理论, 可以保证该方法的局部收敛性:

      $$\left\| {\mathit{\boldsymbol{I}} - \gamma \left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)\mathit{\boldsymbol{G}}\left( {\mathit{\boldsymbol{\hat x}}} \right)\mathit{\boldsymbol{H}}\left( {\mathit{\boldsymbol{\hat x}}} \right)} \right\| < 1$$ (8)

      其中, 海森矩阵为:

      $$\begin{array}{l} \mathit{\boldsymbol{H}}\left( {\mathit{\boldsymbol{\hat x}}} \right) = {\mathit{\boldsymbol{J}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\hat x}}} \right)\mathit{\boldsymbol{PJ}}\left( {\mathit{\boldsymbol{\hat x}}} \right) + \mathop \sum \limits_{i = 1}^n {V_i}\left( {\hat x} \right){p_i}{V_i}\left( {\hat x} \right) = \\ \;\;\;\;\mathit{\boldsymbol{N}}\left( {\mathit{\boldsymbol{\hat x}}} \right) + \mathit{\boldsymbol{S}}\left( {\mathit{\boldsymbol{\hat x}}} \right) > 0 \end{array}$$ (9)

      式中, 由于海森矩阵储存成本高, 计算量比较大, 当处理残差较小、非线性强度较弱的非线性模型时, 可将其近似为$\mathit{\boldsymbol{H}}\left( {\mathit{\boldsymbol{\hat x}}} \right) \approx \mathit{\boldsymbol{N}}\left( {\mathit{\boldsymbol{\hat x}}} \right) = {\mathit{\boldsymbol{J}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\hat x}}} \right)\mathit{\boldsymbol{PJ}}\left( {\mathit{\boldsymbol{\hat x}}} \right)$, $\mathit{\boldsymbol{N}}\left( \mathit{\boldsymbol{x}} \right)$为最小二乘目标函数线性化后系数矩阵。若选取$\lambda \mathit{\boldsymbol{Q}}\left( {\mathit{\boldsymbol{\hat x}}} \right) = {\mathit{\boldsymbol{N}}^{ - 1}}\left( {\mathit{\boldsymbol{\hat x}}} \right)$, 则可构造高斯牛顿法:

      $${\mathit{\boldsymbol{\hat x}}_{k + 1}} = {\mathit{\boldsymbol{\hat x}}_k} + {\mathit{\boldsymbol{N}}^{ - 1}}\left( {\mathit{\boldsymbol{\hat x}}} \right){\mathit{\boldsymbol{J}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\hat x}}} \right)\mathit{\boldsymbol{PV}}\left( {\mathit{\boldsymbol{\hat x}}} \right)$$ (10)

      通常而言, 高斯牛顿法具有线性收敛效率。若残差较大或者非线性强度过高, 高斯牛顿法将会失效, 若伴随不适定问题, 会进一步降低参数估计解的收敛稳定性。牛顿法和高斯牛顿法共同的最大问题是, 当观测模型病态时, 解会出现迭代序列不稳定特征, 甚至无法收敛而解算失效。

    • 非线性测距定位方程最小二乘解的性质是观测向量末端以观测权重为质量的质点系的重心。文献[13]根据非线性最小二乘解性质导出了测距定位方程参数估计的重心迭代法, 其迭代公式为:

      $${\mathit{\boldsymbol{\hat x}}_{k + 1}} = \mathop \sum \limits_{i = 1}^n \frac{{{\mathit{\boldsymbol{S}}_i}\left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)}}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}$$ (11)

      式中, $\boldsymbol{S}_{i}\left(\hat{\boldsymbol{x}}_{k}\right)=p_{i}\left(\boldsymbol{x}_{i}-\boldsymbol{e}_{i}^{k} L_{i}\right) ; \hat{\boldsymbol{x}}_{k+1}$是由解${\mathit{\boldsymbol{\hat x}}_k}$末端以观测权重为质量的质点系的重心, 式(11)也可写为:

      $${\mathit{\boldsymbol{\hat x}}_{k + 1}} = {\mathit{\boldsymbol{\hat x}}_k} + \frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathit{\boldsymbol{J}}{({\mathit{\boldsymbol{\hat x}}_k})^{\rm{T}}}\mathit{\boldsymbol{PV}}\left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)$$ (12)

      式(11)称为重心迭代公式, 而式(12)则表明其实际是一种最速下降法, 即选取$\gamma \left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)\mathit{\boldsymbol{G}}\left( {\mathit{\boldsymbol{\hat x}}} \right) = 1/{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)$。

      重心迭代法不依赖于初始值的精度, 具有计算简单, 无需矩阵取逆、计算海森矩阵的优点。然而, 重心迭代法具有较低的收敛效率, 尤其受到1/tr(P)值的影响。1/tr(P)是重心法求解的保守策略, 易受到观测数据个数n的影响, 即随着观测数目增多, 其线性搜索策略可能越保守。1/tr(P值随观测数据个数n增加而逐渐下降, 并存在趋于零的倾向。例如, 多模全球导航卫星系统(Global Navigation Satellite System, GNSS)融合定位已成为GNSS定位技术发展的重要方向。多模GNSS融合定位中, 单历元采集的可用卫星数能达到几十颗。若采用重心迭代法来处理多模GNSS测距定位方程, 由于1/tr(P值过小, 导致相邻迭代序列步长较短, 收敛效率低且实际应用价值不高。再例如, 在移动测量定位和水下定位领域, 解决三维定位问题时, 采集数据个数需满足n≥4,则1/tr(P)≤0.25。即使给定的1/tr(P能够保证重心迭代法稳定收敛至预期解, 但由于初始精度、模型态性和1/tr(P的影响, 具有较低的收敛速度。为增加重心法的实用性, 必须提高其收敛效率。

    • 1967年, 萨玛斯基将萨玛斯基技巧应用于牛顿法, 综合牛顿法的收敛速度和简化牛顿法计算量小的特点, 将t步牛顿法简化为一步牛顿法, 迭代法具有t+1阶收敛率, 明显提高了牛顿法的收敛效率[15-16]。本文将萨玛斯技巧应用到Barycentre迭代法中, 减少导数计算量, 提高Barycentre迭代法的收敛效率, 称之为Frozen-Barycentre迭代法。该方法采用不动点理论, 在外部迭代基础上加入内部迭代。

      当外部迭代次数k为1时, 固定不动点${\mathit{\boldsymbol{\hat x}}_0}$, 可进行内部迭代, 其迭代序列为:

      $$\left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\hat x}}}_0} = {{\mathit{\boldsymbol{\hat x}}}_0}}\\ {{{\mathit{\boldsymbol{\hat x}}}_{0, 1}} = {{\mathit{\boldsymbol{\hat x}}}_0} + \frac{1}{{{\mathop{\rm tr}\nolimits} (\mathit{\boldsymbol{P}})}}\mathit{\boldsymbol{J}}{{\left( {{{\mathit{\boldsymbol{\hat x}}}_0}} \right)}^{\rm{T}}}\mathit{\boldsymbol{PV}}\left( {{{\mathit{\boldsymbol{\hat x}}}_0}} \right)}\\ {{{\mathit{\boldsymbol{\hat x}}}_{0, 2}} = {{\mathit{\boldsymbol{\hat x}}}_{0, 1}} + \frac{1}{{{\mathop{\rm tr}\nolimits} (\mathit{\boldsymbol{P}})}}\mathit{\boldsymbol{J}}{{\left( {{{\mathit{\boldsymbol{\hat x}}}_0}} \right)}^{\rm{T}}}\mathit{\boldsymbol{PV}}\left( {{{\mathit{\boldsymbol{\hat x}}}_{0, 1}}} \right)}\\ \;\;\;\;\vdots \\ {{{\mathit{\boldsymbol{\hat x}}}_{0, i}} = {{\mathit{\boldsymbol{\hat x}}}_{0, i}} + \frac{1}{{{\mathop{\rm tr}\nolimits} (\mathit{\boldsymbol{P}})}}\mathit{\boldsymbol{J}}{{\left( {{{\mathit{\boldsymbol{\hat x}}}_0}} \right)}^{\rm{T}}}\mathit{\boldsymbol{PV}}\left( {{{\mathit{\boldsymbol{\hat x}}}_{0, i - 1}}} \right)}\\ {{{\mathit{\boldsymbol{\hat x}}}_1} = {{\mathit{\boldsymbol{\hat x}}}_{0, i}}} \end{array}} \right.$$ (13)

      同理, 若$k = 1, 2 \cdots $时, 固定不动点${\mathit{\boldsymbol{\hat x}}_k}$, 其内部迭代序列为:

      $$\left\{ {\begin{array}{*{20}{l}} {{x_{k, 1}} = {x_k} + \frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathit{\boldsymbol{J}}{{({{\mathit{\boldsymbol{\hat x}}}_k})}^{\rm{T}}}\mathit{\boldsymbol{PV}}\left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)}\\ {{x_{k, 2}} = {x_{k, 1}} + \frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathit{\boldsymbol{J}}{{({{\mathit{\boldsymbol{\hat x}}}_k})}^{\rm{T}}}\mathit{\boldsymbol{PV}}\left( {{{\mathit{\boldsymbol{\hat x}}}_{k, 1}}} \right)}\\ {\begin{array}{*{20}{c}} {}&{} \end{array} \vdots }\\ {{x_{k, i}} = {x_{k, i - 1}} + \frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathit{\boldsymbol{J}}{{({{\mathit{\boldsymbol{\hat x}}}_k})}^{\rm{T}}}\mathit{\boldsymbol{PV}}\left( {{{\mathit{\boldsymbol{\hat x}}}_{k, i - 1}}} \right)}\\ {{x_{k + 1}} = {x_{k, i}}} \end{array}} \right.$$ (14)

      由式(13)和式(14)发现, Frozen-Barycentre迭代法是在外部迭代基础上加入内部迭代序列, 通过减少导数计算量来提高重心迭代的收敛效率。具体迭代公式为:

      $$\left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\hat x}}}_{k, i + 1}} = {{\mathit{\boldsymbol{\hat x}}}_{k, i}} + \frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathit{\boldsymbol{J}}{{({{\mathit{\boldsymbol{\hat x}}}_k})}^{\rm{T}}}\mathit{\boldsymbol{PV}}\left( {{{\mathit{\boldsymbol{\hat x}}}_{k, i}}} \right)}\\ {{{\mathit{\boldsymbol{\hat x}}}^k} = {{\mathit{\boldsymbol{\hat x}}}_{k, 0}}, {{\mathit{\boldsymbol{\hat x}}}_{k + 1}} = {{\mathit{\boldsymbol{\hat x}}}_{k, m}}} \end{array}} \right.$$ (15)

      式(15)称为Frozen-Barycentre迭代公式, 其中, k表示外部迭代次数; t表示内部迭代次数, 且t≥1。该方法的优点在于, 内部迭代过程中保持迭代矩阵一致性, 通过降低残差量来提高重心法的收敛效率。这是因为在大多数情况下, 求解非线性系统迭代矩阵$\mathit{\boldsymbol{J}}\left( {\mathit{\boldsymbol{\hat x}}} \right)$的成本比降低残差量$\mathit{\boldsymbol{V}}\left( {\mathit{\boldsymbol{\hat x}}} \right)$的成本要高。当t=1时, Frozen-Barycentre迭代法将退化为Barycentre法。

    • 由式(11)可知, $\mathit{\boldsymbol{\hat x}}$的解析表达式可写为:

      $$\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\hat x}} = \frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathop \sum \limits_{i = 1}^n {p_i}\left( {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{e}}_i}{L_i}} \right) = }\\ {\frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathop \sum \limits_{i = 1}^n {p_i}\left( {{\mathit{\boldsymbol{x}}_i} - \frac{{{\mathit{\boldsymbol{x}}_i} - {{\mathit{\boldsymbol{\hat x}}}_i}}}{{{\mathit{\boldsymbol{x}}_i} - {{\mathit{\boldsymbol{\hat x}}}_i}}}{L_i}} \right)} \end{array}$$ (16)

      结合式(16), 初始值可通过式(17)计算:

      $$\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\hat x}}}_0} = \frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathop \sum \limits_{i = 1}^n {p_i}\left( {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{e}}_i}{L_i}} \right) = }\\ {\frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathop \sum \limits_{i = 1}^n {p_i}\left( {{\mathit{\boldsymbol{x}}_i} - \frac{{{\mathit{\boldsymbol{x}}_i}}}{{{\mathit{\boldsymbol{x}}_i}}}{L_i}} \right)} \end{array}$$ (17)

      其中, $\frac{1}{{{\rm{tr}}\left( \mathit{\boldsymbol{P}} \right)}}\mathop \sum \limits_{i = 1}^n {p_i}\left( {{\mathit{\boldsymbol{x}}_i} - \frac{{{\mathit{\boldsymbol{x}}_i}}}{{{\mathit{\boldsymbol{x}}_i}}}{L_i}} \right)$是已知点相对于坐标原点组成的质点系的重心。

      非线性测距定位方程线性化初值并没有很好的求解策略, 一般采用先验初值或者设初值为零, 先验初值精度相对较高, 但获取较难且成本高。采用式(17)计算初值能够节省先验初值获取成本, 并且相对于零初值, 该初值具有物理几何意义, 可信度高。然而, 水面GNSS浮标和水下应答器点与坐标原点相隔距离较远, 采用式(17)计算的初值精度相对较差。短程水下测距定位方程初始雅克比矩阵奇异或者病态, 采用该初值会导致牛顿型数值方法计算不稳定, 甚至无法收敛; 此外, 也会降低重心迭代法的收敛速度。在GNSS定位中, 用式(17)计算的初值点能够落在控制图形内, 从而加快牛顿型数值方法和重心法的收敛速度, 此时能够利用简单的几步迭代获取一定精度的解, 具有很好的实用价值。

    • 本文采用文献[13]模拟的短程测距定位问题。设${P_1}, {P_2} \ldots {P_8}$为8个已知点, 且已知点到待测点P9(模拟值为(0, 0, 0)m)的观测距离如表 1所示, 距离为等精度观测, 要求根据8个观测距离确定待测点的坐标。迭代终止条件为$\left\| {\mathit{\boldsymbol{h}}\left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)} \right\| \le 1 \times {10^{ - 8}}$。迭代初始值由式(17)计算。

      表 1  短程距离观测控制点坐标和观测距离/m

      Table 1.  Trilateration Network Known Coordinates and Observation Distances/m

      点号 坐标 距离
      X Y Z
      P1 10.000 20.000 6.000 22.180
      P2 8.000 30.000 0.000 30.960
      P3 10.000 -20.000 6.000 23.170
      P4 8.000 -30.000 0.000 30.880
      P5 -10.000 -20.000 -6.000 23.070
      P6 -8.000 -30.000 0.000 30.440
      P7 -10.000 20.000 -6.000 22.700
      P8 -8.000 30.000 0.000 30.720

      该算例中非线性代数方程式(1)共有8个, 其中未知数3个。由非线性代数方程构建的非线性方程组依据泰勒公式线性化后, 由于已知点和未知点近似共面, 方程组系数矩阵具有严重的病态性。由表 1数据初步分析可知, 方程组的病态主要集中在Z方向, X方向和Y方向只存在微弱的病态性。图 1给出了高斯牛顿法的点位迭代序列。由图 1可知, 高斯牛顿法点位迭代序列值所受观测扰动影响较大而无法收敛, 这主要是由于给定的初始值精度差和方程组系数矩阵具有较强的病态性, 导致迭代序列计算极其不稳定以至于无法收敛。本文建议病态方程组数值解算时不采用高斯牛顿法。

      图  1  高斯—牛顿法的点位迭代序列

      Figure 1.  The Sequence of Gauss-Newton Method

      表 2给出了线性化平差、Barycentre和Frozen-Barycentre的解算结果。由表 2可知, 由于方程组病态和线性化初值精度的影响, 短程测距定位方程若直接采用线性化平差估计, 解不稳定且精度差; Frozen-Barycentre和Barycentre法数值收敛解相同, 均提高线性化平差估计解的精度, 这是由于两种方法通过避免矩阵求逆运算, 解决了短程测距定位方程的不适定问题, 确保了解的收敛稳定性; 通过对比Barycentre法的数值收敛解和线性化平差估计解的残差, 表明迭代估计在非线性模型病态求解问题中更加可靠准确。

      表 2  不同方法的解算结果

      Table 2.  Calculation Results of Different Methods

      方法 t 数值收敛解$\mathit{\boldsymbol{\hat x}}$/m 时间/s 外循环k 残差$\left\| {\Delta \mathit{\boldsymbol{\hat x}}} \right\|/{\bf{m}}$
      线性化平差 -0.648 7 0.1237 1.4791 1.6198
      Barycentre -0.214 1 0.1238 0.5493 0.130909 987 0.6024
      Frozen-Barycentre 2 -0.214 1 0.1238 0.5493 0.119818 491 0.6024
      3 -0.214 1 0.1238 0.5493 0.097318 326 0.6024
      4 -0.214 1 0.1238 0.5493 0.085686 243 0.6024
      5 -0.214 1 0.1238 0.5493 0.078891 193 0.6024
      6 -0.214 1 0.1238 0.5493 0.073928 160 0.6024
      7 -0.214 1 0.1238 0.5493 0.073005 137 0.6024
      8 -0.214 1 0.1238 0.5493 0.067277 119 0.6024
      9 -0.214 1 0.1238 0.5493 0.065978 105 0.6024
      10 -0.214 1 0.1238 0.5493 0.064174 94 0.6024

      Frozen-Barycentre法加入内部迭代序列并不会影响算法的最终收敛结果; 该方法外部迭代终止步k虽然低于Barycentre法, 但系统由于内部迭代t的选择并没有节省整体循环次数, 甚至有所提高, 这是由于Frozen-Barycentre法是在Barycentre法的基础上加入萨玛斯基技巧, 通过内循环和外循环的转换来减少导数计算量, 节省算法的运算时间。由表 2能够看出, Frozen-Barycentre迭代法循环消耗时间优于Barycentre法, 表明Frozen-Barycentre迭代法具有更好的收敛效率。经过模拟数据研究, 内部循环迭代次数设置为10次。

      为比较Frozen-Barycentre法和Barycentre法的点位迭代序列, 本文设置截断次数为200次。图 2给出了两种方法的点位迭代序列。由图 2可知, Frozen-Barycentre法继承了Barycentre迭代法收敛序列稳定的优点, 明显提高了Barycentre法的收敛效率。与图 1相比可知, Frozen-Barycentre法和Barycentre法通过避免矩阵取逆, 解决了方程组的病态问题, 能够稳定收敛至最小二乘解。此外, Frozen-Barycentre法和Barycentre法具有对初值精度依赖性低的优点。

      图  2  Barycentre和Frozen-Barycentre法的点位迭代序列

      Figure 2.  The Sequences of Barycentre and Frozen-Barycentre Methods

    • 为进一步研究Barycentre法和Frozen-Barycentre法的实用性, 本文采用南海实测水下定位数据, 应答器位于水下2 000 m左右, 一个载有USBL的测量船围绕海底应答器航行, 共采集了724个观测数据, 测量分布图如图 3所示。在该数据中选取15个相邻数据进行验证, 数据坐标如图 4所示。假设观测距离为等精度观测, 要求根据15个已知点坐标和观测距离求待测点的坐标。以全部724个控制点观测数据综合定位结果作为真值。迭代终止条件为$\left\| {\mathit{\boldsymbol{h}}\left( {{{\mathit{\boldsymbol{\hat x}}}_k}} \right)} \right\| \le 1 \times {10^{ - 8}}$。数据采集时, 测量船在水平面航行, Z方向上具有较强的病态性; 此外, 测量船航行过程中, 采集数据相隔时间较短, X方向和Y方向也会产生微弱的病态性。

      图  3  测量船围绕应答器实测航线

      Figure 3.  The Ship's Measured Route Around the Transponder

      图  4  15个相邻数据的已知点坐标

      Figure 4.  Known Point Coordinates for 15 Adjacent Data

      图 4和已知点数据初步分析可知, 所构建非线性方程组依据泰勒公式转为线性方程组后, 由于未知点与已知点近似共面, 导致方程组系数矩阵具有严重的病态, 并且病态性主要集中在Z方向。采用式(17)计算初始值。该案例中高斯牛顿法由于初始值精度和方程组病态的影响而无法收敛, 解算失效。

      表 3给出了线性化平差、Barycentre法和Frozen-Barycentre法的解算结果。由于水下短程测距定位方程系数矩阵具有较强的病态性和线性化初值精度的影响, 直接采用线性化平差估计, 其估计解与真解相差较大, 解精度低且不稳定。Frozen-Barycentre法和Barycentre法的解算结果一致, 均明显提高了线性化平差估计解的精度。

      表 3  不同算法的水下定位解算结果

      Table 3.  Results of Convergence of Barycentre and Frozen-Barycentre Algorithms

      方法 t 数值收敛解/m 时间/s k 残差$\left\| {\Delta \mathit{\boldsymbol{\hat x}}} \right\|/{\rm{m}}$
      线性化平差 2 553993.560 3 -145 334.776 5 371795 350.815 7 371793 904.189
      Barycentre 2 438949.021 9 492011.361 4 2142.926 3 74.9428 25 032 135.102
      Frozen-Barycentre 2 2 438949.021 9 492011.361 4 2142.926 3 20.2227 12 518 135.102
      3 2 438949.021 9 492011.361 4 2142.926 3 10.3040 8 347 135.102
      4 2 438949.021 9 492011.361 4 2142.926 3 7.2169 6 262 135.102
      5 2 438949.021 9 492011.361 4 2142.926 3 5.0740 5 010 135.102
      6 2 438949.021 9 492011.361 4 2142.926 3 4.3751 4 176 135.102
      7 2 438949.021 9 492011.361 4 2142.926 3 3.5661 3 580 135.102
      8 2 438949.021 9 492011.361 4 2142.926 3 3.2884 3 133 135.102
      9 2 438949.021 9 492011.361 4 2142.926 3 2.9414 2 786 135.102
      10 2 438949.021 9 492011.361 4 2142.926 3 2.7825 2 508 135.102

      Barycentre法通过避免矩阵取逆运算, 能够解决水下短程测距定位方程的不适定问题, 但收敛效率极低, 尤其是运算时间。Frozen-Barycentre迭代法在重心法基础上加入萨玛斯基技巧, 实现了内迭代与外迭代的转换, 通过减少导数计算量, 节省了Barycentre法的运算时间, 具有更好的收敛效率, 实用价值更高。当然, Frozen-Barycentre迭代法的整体循环次数并未有明显的改善, 甚至有所提升, 比如内部迭代次数为10次时, Frozen-Barycentre整体循环次数为25 080次, 高于25 032次, 但是计算机相对算法整体循环次数更加注重运算时间。实测数据验证, 考虑内迭代和外迭代成本, 内迭代循环次数一般选择为10次。

      为说明Frozen-Barycentre法和Barycentre法的点位迭代序列变化情况, 图 5给出了两种方法在截断次数为9 000次的点位迭代序列(迭代次数9 000次后, Barycentre法的点位迭代序列比较平稳)。由图 5可知, Frozen-Barycentre法和Barycentre在X方向和Y方向点位迭代序列变化较为平缓; 在Z方向上, Frozen-Barycentre法和Barycentre法点位迭代序列整体变化较为平缓, 可是分别在迭代341次和3 311次都产生了较大的扰动, 这是由于短程水下测距定位方程病态性主要集中于Z方向, 导致迭代序列产生不稳定特征; Frozen-Barycentre法通过内部迭代虽加快了外部迭代进程, 但与Barycentre法点位迭代序列变化趋势几乎一致, 表明该方法在继承重心法的迭代稳定的优点的同时提高了Barycentre法的收敛效率。

      图  5  Barycentre和Frozen-Barycentre法的点位迭代序列

      Figure 5.  The Sequences of Barycentre and Frozen-Barycentre Methods

    • 高斯牛顿法处理短程测距定位问题时, 由于受到初始值精度和方程组病态的影响而无法收敛。Barycentre法通过避免矩阵求逆运算, 能够稳定收敛至最小二乘解, 但Barycentre法受到步长因子1/tr(P)值的影响, 收敛效率较低。本文将萨玛斯基应用于Barycentre法, 提出了Frozen-Barycentre迭代法, 该方法继承了Barycentre法不依赖于初始值的精度, 具有计算简单、无需矩阵求逆、计算海森矩阵的优点, 通过内迭代和外迭代的转换减少导数计算量来节省算法的运算时间, 提高Barycentre法的收敛速度。算例结果表明, 线性化平差解精度低且误差大; 高斯-牛顿法由于初始值精度和方程组系数矩阵病态而无法收敛; 新方法明显提高了线性化平差估计解的解算精度, 且收敛效率应优于Barycentre方法, 具有更好的实用价值。

参考文献 (16)

目录

    /

    返回文章
    返回