留言板

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

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

三维坐标转换的高斯-赫尔默特模型及其抗差解法

刘超 王彬 赵兴旺 余学祥

刘超, 王彬, 赵兴旺, 余学祥. 三维坐标转换的高斯-赫尔默特模型及其抗差解法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348
引用本文: 刘超, 王彬, 赵兴旺, 余学祥. 三维坐标转换的高斯-赫尔默特模型及其抗差解法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348
LIU Chao, WANG Bin, ZHAO Xingwang, YU Xuexiang. Three-Dimensional Coordinate Transformation Model and Its Robust Estimation Method Under Gauss-Helmert Model[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348
Citation: LIU Chao, WANG Bin, ZHAO Xingwang, YU Xuexiang. Three-Dimensional Coordinate Transformation Model and Its Robust Estimation Method Under Gauss-Helmert Model[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348

三维坐标转换的高斯-赫尔默特模型及其抗差解法

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

国家自然科学基金 41404004

国家自然科学基金 41474026

安徽省博士后基金 2015B044

安徽理工大学科研启动基金 11152

详细信息
    作者简介:

    刘超, 博士, 副教授, 主要从事变形数据分析理论与方法研究。chaoliu0202@gmail.com

  • 中图分类号: P207

Three-Dimensional Coordinate Transformation Model and Its Robust Estimation Method Under Gauss-Helmert Model

Funds: 

The National Natural Science Foundation of China 41404004

The National Natural Science Foundation of China 41474026

the Postdoctoral Science Foundation of Anhui Province of China 2015B044

the Anhui University of Science and Technology Foundation 11152

More Information
    Author Bio:

    LIU Chao, PhD, associate professor, specializes in deformation monitoring data analysis. E-mail:chaoliu0202@gmail.com

  • 摘要: 对三维坐标转换的高斯-赫尔默特(Gauss-Helmert,GH)模型,采用牛顿-高斯(Newton-Gauss)迭代算法构建了该模型的拉格朗日目标函数,推导了其解算方法,并给出了具体的计算步骤。在此基础上,考虑到可能出现的粗差对观测空间与结构空间的综合影响,基于标准化残差构造权因子函数,推导了该模型的抗差解法。仿真实验结果表明,GH模型用于三维坐标转换时不受旋转角度大小和其他附加条件限制,解算结果与现有算法一致,且估计参数的维数大大降低,计算效率有一定程度的提高;所提出的抗差解法效果良好,与现有基于整体最小二乘的三维坐标转换的抗差解法相比,表现出了更好的稳健性。
  • 图  1  方案5和6估值与真值的差值序列(3个粗差)

    Figure  1.  Difference Values Between Truth Values and Estimated Results of Schemes 5 and 6 with Three Outliers

    表  1  存在1个粗差时各方案的均方根误差对比

    Table  1.   Root Mean Square Errors of Estimated Parameters for Different Schemes with One Outlier

    方案X0 /mY0 /mZ0 /mμ/(10-5 rad)α1 /(10-5 rad)α2 /(10-5 rad)α3 /(10-5 rad)
    方案119.302 110.579 615.580 12.6001.5101.1601.750
    方案23.922 12.666 33.613 20.6320.4050.3510.366
    方案33.92 212.666 33.613 20.6320.4050.3510.366
    方案46.470 33.981 95.938 81.0100.6440.5570.570
    方案54.502 93.163 64.381 60.6960.4530.4100.452
    方案64.194 92.907 24.106 60.6660.4350.3820.407
    下载: 导出CSV

    表  2  存在3个粗差时各方案的均方根误差对比

    Table  2.   Root Mean Square Errors of Estimated Parameters for Different Schemes with Three Outliers

    方案X0 /mY0 /mZ0 /mμ/(10-5 rad)α1 /(10-5 rad)α2 /(10-5 rad)α3 /(10-5 rad)
    方案110.307 912.340 718.474 01.4601.3901.7301.130
    方案23.752 92.956 42.844 90.5070.2890.3110.417
    方案33.752 92.956 42.844 90.5070.2890.3110.417
    方案48.797 86.880 56.709 91.2300.6640.7700.952
    方案57.068 05.744 25.810 60.9780.5920.6280.798
    方案64.695 83.896 13.931 20.6710.4040.4130.518
    下载: 导出CSV

    表  3  存在5个粗差时各方案的均方根误差对比

    Table  3.   Root Mean Square Errors of Estimated Parameters for Different Schemes with Five Outliers

    方案X0 /mY0 /mZ0 /mμ/(10-5 rad)α1 /(10-5 rad)α2 /(10-5 rad)α3 /(10-5 rad)
    方案112.748 611.219 615.764 92.2501.7901.1501.400
    方案23.414 03.380 14.576 30.6100.4690.3800.424
    方案33.414 03.380 14.576 30.6100.4690.3800.424
    方案410.173 610.600 614.184 51.9601.4801.1201.320
    方案58.900 88.467 611.919 81.5801.2900.9561.110
    方案66.064 66.039 37.999 31.0900.8180.6390.757
    下载: 导出CSV

    表  4  方案5和方案6计算检查点的均方根误差(1个粗差)/m

    Table  4.   Root Mean Square Errors of Check Points for Schemes 5 and 6 with One Outlier/m

    方向方案点位1点位2点位3点位4点位5点位6点位7
    X方案50.079 20.081 30.072 90.077 60.081 50.080 30.080 3
    方案60.078 40.080 30.072 40.077 10.080 80.079 00.078 7
    Y方案50.096 20.040 70.076 90.039 80.065 60.068 20.041 0
    方案60.095 90.039 10.077 00.038 60.065 30.068 20.040 1
    Z方案50.084 80.052 20.067 50.059 80.077 60.082 80.057 3
    方案60.084 70.052 20.066 80.058 30.076 40.081 90.056 1
    下载: 导出CSV

    表  5  方案5和方案6计算检查点的均方根误差(3个粗差)/m

    Table  5.   Root Mean Square Errors of Check Points for Schemes 5 and 6 with Three Outliers/m

    方向方案点位1点位2点位3点位4点位5点位6点位7
    X方案50.078 50.081 10.073 90.082 20.082 70.080 30.082 9
    方案60.077 80.078 50.071 70.078 20.080 40.078 50.078 0
    Y方案50.098 70.048 10.082 70.044 60.065 10.073 00.044 9
    方案60.094 60.041 90.078 40.038 70.063 90.070 70.039 2
    Z方案50.089 30.056 90.068 30.061 80.076 80.083 40.059 3
    方案60.087 10.052 90.066 20.057 60.076 20.082 10.057 2
    下载: 导出CSV

    表  6  方案5和方案6计算检查点的均方根误差(5个粗差)/m

    Table  6.   Root Mean Square Errors of Check Points for Schemes 5 and 6 with Five Outliers/m

    方向方案点位1点位2点位3点位4点位5点位6点位7
    X方案50.084 40.09 220.080 60.089 30.087 30.090 70.097 6
    方案60.079 70.083 50.075 00.082 10.081 30.081 60.085 5
    Y方案50.104 00.057 70.087 70.052 20.072 80.079 40.057 6
    方案60.097 90.047 50.083 10.041 90.069 50.072 10.046 8
    Z方案50.099 50.066 10.078 10.073 40.085 40.087 90.068 1
    方案60.091 90.060 10.070 80.062 10.078 90.083 50.060 0
    下载: 导出CSV
  • [1] Grafarend E W. Nonlinear Analysis of the Three-Dimensional Datum Transformation[J].Journal of Geodesy, 2003, 77(1):66-76 doi:  10.1007-s00190-002-0299-9/
    [2] 姚宜斌, 黄承猛, 李程春, 等.一种适用于大角度三维坐标转换参数求解算法[J].武汉大学学报·信息科学版, 2012, 37(3):253-256 http://ch.whu.edu.cn/CN/Y2012/V37/I3/253

    Yao Yibin, Huang Chengmeng, Li Chengchun, et al.A New Algorithm for Solution of Transformation Parameters of Big Rotation Angle's 3D Coordinate[J].Geomatics and Information Science of Wuhan University, 2012, 37(3):253-256 http://ch.whu.edu.cn/CN/Y2012/V37/I3/253
    [3] 杨元喜, 徐天河.不同坐标系综合变换法[J].武汉大学学报·信息科学版, 2001, 26(6):509-513 http://ch.whu.edu.cn/CN/abstract/abstract5223.shtml

    Yang Yuanxi, Xu Tianhe.The Combined Method of Datum Transformation Between Different Coordinate Systems[J].Geomatics and Information Science of Wuhan University, 2001, 26(6):509-513 http://ch.whu.edu.cn/CN/abstract/abstract5223.shtml
    [4] 刘大杰, 施一民, 过静珺.全球定位系统(GPS)的原理与数据处理[M].上海:同济大学出版社, 1996

    Liu Dajie, Shi Yimin, Guo Jingjun.The Theory and Data Processing of Global Positioning System (GPS)[M].Shanghai:Tongji University Press, 1996
    [5] 曾文宪, 陶本藻.三维坐标转换的非线性模型[J].武汉大学学报·信息科学版, 2003, 28(5):566-568 http://ch.whu.edu.cn/CN/abstract/abstract4861.shtml

    Zeng Wenxian, Tao Benzao.Non-linear Adjustment Model of Three-Dimensional Coordinate Transformation[J].Geomatics and Information Science of Wuhan University, 2003, 28(5):566-568 http://ch.whu.edu.cn/CN/abstract/abstract4861.shtml
    [6] 陆珏, 陈义, 郑波.总体最小二乘方法在三维坐标转换中的应用[J].大地测量学与地球动力学, 2008, 28(5):77-81 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200805016

    Lu Jue, Chen Yi, Zheng Bo. Total Least Squares to Three-Dimensional Datum Transformation[J].Journal of Geodesy and Geodynamics, 2008, 28(5):77-81 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200805016
    [7] Golub H G, Vanloan F C.An Analysis of the Total Least Squares Problem[J].SIAM Journal on Numerical Analysis, 1980, 17(6):883-893 doi:  10.1137/0717073
    [8] Schaffrin B, Felus Y A.On the Multivariate Total Least-Squares Approach to Empirical Coordinate Transformations[J].Journal of Geodesy, 2008, 82(6):373-383 doi:  10.1007/s00190-007-0186-5
    [9] Felus F, Burtch R.On Symmetrical Three-Dimensional Datum Conversion[J]. GPS Solutions, 2009, 13(1):65-74 doi:  10.1007/s10291-008-0100-5
    [10] Shen Yunzhong, Li Bofeng, Chen Yi. An Iterative Solution of Weighted Total Least-Squares Adjustment[J]. Journal of Geodesy, 2011, 85(4):229-238 doi:  10.1007/s00190-010-0431-1
    [11] Xu Peiliang, Liu Jingnan, Shi Chuang.Total Least Squares Adjustment in Partial Errors-in-Variables Models:Algorithm and Statistical Analysis[J].Journal of Geodesy, 2012, 86(8):661-675 doi:  10.1007/s00190-012-0552-9
    [12] Xu Peiliang, Liu Jingnan.Variance Components in Errors-in-Variables Models:Estimability, Stability and Bias Analysis[J].Journal of Geodesy, 2014, 88(8):719-734 doi:  10.1007/s00190-014-0717-9
    [13] 魏二虎, 殷志祥, 李广文, 等.虚拟观测值法在三维坐标转换中的应用研究[J].武汉大学学报·信息科学版, 2014, 39(2):152-156 http://ch.whu.edu.cn/CN/abstract/abstract2918.shtml

    Wei Erhu, Yin Zhixiang, Li Guangwen, et al.On 3D Coordinate Transformation with Virtual Observation Method[J].Geomatics and Information Science of Wuhan University, 2014, 39(2):152-156 http://ch.whu.edu.cn/CN/abstract/abstract2918.shtml
    [14] 方兴, 曾文宪, 刘经南, 等.三维坐标转换的通用整体最小二乘算法[J].测绘学报, 2014, 43(11):1139-1143 http://www.cqvip.com/QK/90069X/201411/663183594.html

    Fang Xing, Zeng Wenxian, Liu Jingnan, et al.A General Total Least Squares Algorithm for Three-Dimensional Coordinate Transformations[J].Acta Geodaetica et Cartographica Sinica, 2014, 43(11):1139-1143 http://www.cqvip.com/QK/90069X/201411/663183594.html
    [15] Neitzel F.Generalization of Total Least-Squares on Example of Unweighted and Weighted 2D Similarity Transformation[J].Journal of Geodesy, 2010, 84(12):751-762 doi:  10.1007/s00190-010-0408-0
    [16] 方兴, 曾文宪, 刘经南, 等.基于非线性高斯-赫尔默特模型的混合整体最小二乘估计[J].测绘学报, 2016, 45(3):291-296 http://d.old.wanfangdata.com.cn/Periodical/chxb201603007

    Fang Xing, Zeng Wenxian, Liu Jingnan, et al.Mixed LS-TLS Estimation Based on Nonlinear Gauss-Helmert Model[J].Acta Geodaetica et Cartographi-ca Sinica, 2016, 45(3):291-296 http://d.old.wanfangdata.com.cn/Periodical/chxb201603007
    [17] Chang Guobin.On Least-Squares Solution to 3D Similarity Transformation Problem Under Gauss-Helmert Model[J].Journal of Geodesy, 2015, 89(6):573-576 doi:  10.1007/s00190-015-0799-z
    [18] 欧吉坤.一种检测粗差的新方法——拟准测定法[J].科学通报, 1999, 44(16):1777-1781 doi:  10.3321/j.issn:0023-074X.1999.16.021

    Ou Jikun.A New Method to Detect Gross Errors:Quasi-Accurate Detection Method[J].Chinese Science Bulletin, 1999, 44(16):1777-1781 doi:  10.3321/j.issn:0023-074X.1999.16.021
    [19] 欧吉坤.粗差的拟准检定法(QUAD法)[J].测绘学报, 1999, 28(1):15-20 doi:  10.3321/j.issn:1001-1595.1999.01.004

    Ou Jikun. Quasi-Accurate Detection of Gross Errors(QUAD)[J].Acta Geodaetica et Cartographica Sinica, 1999, 28(1):15-20 doi:  10.3321/j.issn:1001-1595.1999.01.004
    [20] Yang Yuanxi.Robust Estimation of Geodetic Datum Transformation[J].Journal of Geodesy, 1999, 73(5):268-274 doi:  10.1007/s001900050243
    [21] 周江文, 黄幼才, 杨元喜, 等.抗差估计最小二乘法[M].武汉:华中理工大学出版社, 1997

    Zhou Jiangwen, Huang Youcai, Yang Yuanxi, et al.Robust Least Squares Method[M].Wuhan:Huazhong University of Science and Technology Press, 1997
    [22] 陈义, 陆珏.以三维坐标转换为例解算稳健总体最小二乘方法[J].测绘学报, 2012, 41(5):715-722 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK201205237128

    Chen Yi, Lu Jue. Peforming 3D Similarity Transformation by Robust Total Least Squares[J].Acta Geodaetica et Cartographica Sinica, 2012, 41(5):715-722 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK201205237128
    [23] Lu Jue, Chen Yi, Li Bofeng, et al.Robust Total Least Squares with Reweighting Iteration for Three-Dimensional Similarity Transformation[J].Survey Review, 2014, 46(334):28-36 doi:  10.1179/1752270613Y.0000000050
    [24] Mahboub V, Amirisimkooei A R, Sharifi M A.Ite-ratively Reweighted Total Least Squares:A Robust Estimation in Errors-in-Variables Models[J].Survey Review, 2013, 45(329):92-99 doi:  10.1179/1752270612Y.0000000017
    [25] 龚循强, 李志林.稳健加权总体最小二乘法[J].测绘学报, 2014, 43(9):888-894 http://d.old.wanfangdata.com.cn/Periodical/kckxjs201301005

    Gong Xunqiang, Li Zhilin.A Robust Weighted Total Least Squares Method[J].Acta Geodaetica et Cartographica Sinica, 2014, 43(9):888-894 http://d.old.wanfangdata.com.cn/Periodical/kckxjs201301005
    [26] 王彬, 李建成, 高井祥, 等.抗差加权整体最小二乘模型的牛顿-高斯算法[J].测绘学报, 2015, 44(6):602-608 http://d.old.wanfangdata.com.cn/Periodical/chxb201506002

    Wang Bin, Li Jiancheng, Gao Jingxiang, et al.Newton-Gauss Algorithm of Robust Weighted Total Least Squares Model[J].Acta Geodaetica et Cartographica Sinica, 2015, 44(6):602-608 http://d.old.wanfangdata.com.cn/Periodical/chxb201506002
    [27] Yang Yuanxi. Robust Estimation for Correlated Observations Based on Bifactor Equivalent Weights[J].Journal of Geodesy, 2002, 76(6):353-358 doi:  10.1007/s00190-002-0256-7
  • [1] 谢建, 龙四春, 周璀.  不等式约束PEIV模型的经典最小二乘方法 . 武汉大学学报 ● 信息科学版, 2021, 46(9): 1291-1297. doi: 10.13203/j.whugis20190196
    [2] 曾文宪, 刘泽邦, 方兴, 李玉兵.  通用EIV平差模型的线性化估计算法 . 武汉大学学报 ● 信息科学版, 2021, 46(9): 1284-1290. doi: 10.13203/j.whugis20200243
    [3] 聂建亮, 刘晓云, 田婕, 李秀明, 赵大江, 黄功文, 张海平.  利用自适应水准网动态平差建立山东垂直运动速度模型 . 武汉大学学报 ● 信息科学版, 2020, 45(4): 620-625. doi: 10.13203/j.whugis20180296
    [4] 王乐洋, 赵英文.  非线性平差精度评定的自适应蒙特卡罗法 . 武汉大学学报 ● 信息科学版, 2019, 44(2): 206-213, 220. doi: 10.13203/j.whugis20170064
    [5] 刘春阳, 王坚, 王彬, 刘超, 刘纪平.  基于中位参数法相关观测的抗差加权整体最小二乘算法 . 武汉大学学报 ● 信息科学版, 2019, 44(3): 378-384. doi: 10.13203/j.whugis20160516
    [6] 薛树强, 杨元喜.  抗差高斯-雅克比组合平差法 . 武汉大学学报 ● 信息科学版, 2015, 40(7): 932-937. doi: 10.13203/j.whugis20130247
    [7] 姚宜斌, 黄书华, 张良, 胡羽丰, 李国平.  求解三维坐标转换参数的整体最小二乘新方法 ! . 武汉大学学报 ● 信息科学版, 2015, 40(7): 853-857. doi: 10.13203/j.whugis20130302
    [8] 姚宜斌, 黄书华, 陈家君.  求解自回归模型参数的整体最小二乘新方法 . 武汉大学学报 ● 信息科学版, 2014, 39(12): 1463-1466.
    [9] 葛旭明, 伍吉仓.  三维基准转换的约束加权混合整体最小二乘的迭代解法 . 武汉大学学报 ● 信息科学版, 2012, 37(2): 178-182.
    [10] 姚宜斌, 黄承猛, 李程春, 孔建.  一种适用于大角度的三维坐标转换参数求解算法 . 武汉大学学报 ● 信息科学版, 2012, 37(3): 253-256.
    [11] 蒋庆仙, 王成宾, 马小辉, 陈晓璧.  利用中位数进行光纤陀螺信号抗差估计 . 武汉大学学报 ● 信息科学版, 2011, 36(6): 656-659.
    [12] 戴吾蛟, 丁晓利, 朱建军.  基于载噪比及参数先验信息的抗差模型 . 武汉大学学报 ● 信息科学版, 2008, 33(8): 834-837.
    [13] 姚吉利.  三维坐标转换的静态滤波模型 . 武汉大学学报 ● 信息科学版, 2005, 30(9): 825-828.
    [14] 杨元喜, 高为广.  基于多传感器观测信息抗差估计的自适应融合导航 . 武汉大学学报 ● 信息科学版, 2004, 29(10): 885-888.
    [15] 曾文宪, 陶本藻.  三维坐标转换的非线性模型 . 武汉大学学报 ● 信息科学版, 2003, 28(5): 566-568.
    [16] 余学祥, 吕伟才.  基于标准化残差的相关观测抗差估计模型 . 武汉大学学报 ● 信息科学版, 1999, 24(1): 75-78.
    [17] 杨世清, 余学祥, 吕伟才.  粗差估值型抗差估计 . 武汉大学学报 ● 信息科学版, 1998, 23(1): 10-13.
    [18] 黄幼才, 黄坚.  人口函数的数据处理方法及其抗差估计的应用 . 武汉大学学报 ● 信息科学版, 1993, 18(3): 95-102.
    [19] 张方仁.  最小二乘估计残差的统计性质和应用 . 武汉大学学报 ● 信息科学版, 1990, 15(2): 91-99.
    [20] 王任享.  选权迭代定位粗差时权函数参数之功能 . 武汉大学学报 ● 信息科学版, 1988, 13(4): 42-50.
  • 加载中
图(1) / 表(6)
计量
  • 文章访问数:  2690
  • HTML全文浏览量:  135
  • PDF下载量:  427
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-12-28
  • 刊出日期:  2018-09-05

三维坐标转换的高斯-赫尔默特模型及其抗差解法

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

    国家自然科学基金 41404004

    国家自然科学基金 41474026

    安徽省博士后基金 2015B044

    安徽理工大学科研启动基金 11152

    作者简介:

    刘超, 博士, 副教授, 主要从事变形数据分析理论与方法研究。chaoliu0202@gmail.com

  • 中图分类号: P207

摘要: 对三维坐标转换的高斯-赫尔默特(Gauss-Helmert,GH)模型,采用牛顿-高斯(Newton-Gauss)迭代算法构建了该模型的拉格朗日目标函数,推导了其解算方法,并给出了具体的计算步骤。在此基础上,考虑到可能出现的粗差对观测空间与结构空间的综合影响,基于标准化残差构造权因子函数,推导了该模型的抗差解法。仿真实验结果表明,GH模型用于三维坐标转换时不受旋转角度大小和其他附加条件限制,解算结果与现有算法一致,且估计参数的维数大大降低,计算效率有一定程度的提高;所提出的抗差解法效果良好,与现有基于整体最小二乘的三维坐标转换的抗差解法相比,表现出了更好的稳健性。

English Abstract

刘超, 王彬, 赵兴旺, 余学祥. 三维坐标转换的高斯-赫尔默特模型及其抗差解法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348
引用本文: 刘超, 王彬, 赵兴旺, 余学祥. 三维坐标转换的高斯-赫尔默特模型及其抗差解法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348
LIU Chao, WANG Bin, ZHAO Xingwang, YU Xuexiang. Three-Dimensional Coordinate Transformation Model and Its Robust Estimation Method Under Gauss-Helmert Model[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348
Citation: LIU Chao, WANG Bin, ZHAO Xingwang, YU Xuexiang. Three-Dimensional Coordinate Transformation Model and Its Robust Estimation Method Under Gauss-Helmert Model[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348
  • 三维坐标转换是测绘科学领域的一个重要问题[1-3],经典方法采用诸如布尔萨-沃尔夫模型(Bursa-Wolf)、莫洛坚斯基模型(Molodensky)和武测模型等[4],建立高斯-马尔科夫模型(Gauss-Markov,GM),基于最小二乘(least squares,LS)准则进行参数估计[5]。虽然采用GM模型进行LS估计可以对参数进行有效估计,但考虑到三维坐标转换本质上是一个拟合问题,其系数矩阵同样存在随机误差,经典LS方法难以获得三维坐标转换参数的最优估值[6]。误差变量模型(errors-in-variables,EIV)及整体最小二乘估计方法(total least squares,TLS)的提出为解决这一问题提供了一条有效途径[7-10]

    引入TLS进行三维坐标转换,采用Bursa-Wolf模型,由于比例因子与旋转矩阵的非线性交互影响,以及系数矩阵中固定列的存在,难以获得其标准EIV模型。通过建立部分误差变量模型(partial errors-in-variables,PEIV)或虚拟观测值模型,可以解决上述问题[11-14]。高斯-赫尔默特(Gauss-Helmert,GH)模型是EIV模型的一般性表达,其对于函数关系的处理不需要强加限制条件,采用GH模型建立Bursa-Wolf模型,理论上同样可以解决上述问题[15-17]。文献[17]给出的GH模型简洁,本文基于牛顿-高斯(Newton-Gauss)迭代算法,推导了其解算方法。

    同时,考虑到观测向量和系数矩阵中可能存在粗差,有必要进一步研究EIV模型中粗差的处理方法。基于LS估计的粗差处理方法总体上可分为两类:(1)将粗差归于函数模型,研究粗差探测方法[18-19];(2)将粗差归于随机模型,研究抗差估计方法[20-21]。将相关研究扩展至TLS领域,学者们也进行了一些有益的探讨,尤其是关于TLS的抗差估计方法——稳健加权整体最小二乘估计(robust weighted total least squares, RWTLS)[22-25]。文献[22-25]主要采用残差进行权因子函数的构造,无法顾及结构空间的影响。文献[26]采用标准化残差构造权因子函数,提出的RWTLS算法实现了对观测空间和结构空间同时抗差,取得了更为理想的效果。本文在推导文献[17]给出的GH模型的基本解法的基础上,借鉴文献[26]的思路,推导该模型的残差协因数阵的表达式,获得标准化残差,并构造权因子函数,实现三维坐标转换参数的抗差估计。

    • 三维坐标转换的GH模型可以表示为[17]

      $$ \left[ {\begin{array}{*{20}{c}} X\\ Y\\ Z \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{e_X}}\\ {{e_Y}}\\ {{e_Z}} \end{array}} \right] = \mu \mathit{\boldsymbol{M}}\left[ {\begin{array}{*{20}{c}} {x - {e_x}}\\ {y - {e_y}}\\ {z - {e_z}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{X_0}}\\ {{Y_0}}\\ {{Z_0}} \end{array}} \right] $$ (1)

      式中,[X  Y  Z]T和[eX  eY  eZ]T分别为目标坐标系下的观测向量和其随机误差;[x  y  z]T和[ex  ey  ez]T分别为原坐标系下的观测向量和其随机误差;[X0  Y0  Z0]T为平移参数; μ为尺度参数;M为旋转矩阵,且

      $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}} = {\mathit{\boldsymbol{M}}_3}{\mathit{\boldsymbol{M}}_2}{\mathit{\boldsymbol{M}}_1} = \left[ {\begin{array}{*{20}{c}} {\cos {\alpha _3}}&{\sin {\alpha _3}}&0\\ { - \sin {\alpha _3}}&{\cos {\alpha _3}}&0\\ 0&0&1 \end{array}} \right] \cdot }\\ {\left[ {\begin{array}{*{20}{c}} {\cos {\alpha _2}}&0&{ - \sin {\alpha _2}}\\ 0&1&0\\ {\sin {\alpha _2}}&0&{\cos {\alpha _2}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&{\cos {\alpha _1}}&{\sin {\alpha _1}}\\ 0&{ - \sin {\alpha _1}}&{\cos {\alpha _1}} \end{array}} \right]} \end{array} $$

      其中,α1α2α3为旋转参数。

      基于Newton-Gauss法,采用非线性最小二乘法,对待估参数(X0, Y0, Z0, μ, α1, α2, α3)进行迭代计算[10]。设第i次迭代计算后获得的参数估值为(X0i, Y0i, Z0i, μi, α1i, α2i, α3i),预测残差为(exi, eyi, ezi)(上标i为迭代次数)。对式(1)等式右端在((X0i, Y0i, Z0i, μi, α1i, α2i, α3i), (exi, eyi, ezi))处进行泰勒级数展开,并保留一阶项:

      $$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} X\\ Y\\ Z \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{e_X}}\\ {{e_Y}}\\ {{e_Z}} \end{array}} \right] = {\mu ^i}{\mathit{\boldsymbol{M}}^i}\left[ {\begin{array}{*{20}{c}} {x - e_x^i}\\ {y - e_y^i}\\ {z - e_z^i} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {X_0^i}\\ {Y_0^i}\\ {Z_0^i} \end{array}} \right] + }\\ {{\mathit{\boldsymbol{M}}^i}\left[ {\begin{array}{*{20}{c}} {x - e_x^i}\\ {y - e_y^i}\\ {z - e_z^i} \end{array}} \right]{\rm{d}}\mu + {\mu ^i}\left( {\frac{{\partial \mathit{\boldsymbol{M}}}}{{\partial {\alpha _1}}}{\rm{d}}{\alpha _1} + \frac{{\partial \mathit{\boldsymbol{M}}}}{{\partial {\alpha _2}}}{\rm{d}}{\alpha _2} + } \right.}\\ {\left. {\frac{{\partial \mathit{\boldsymbol{M}}}}{{\partial {\alpha _3}}}{\rm{d}}{\alpha _3}} \right)\left[ {\begin{array}{*{20}{c}} {x - e_x^i}\\ {y - e_y^i}\\ {z - e_z^i} \end{array}} \right] - {\mu ^i}{\mathit{\boldsymbol{M}}^i}\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{e_x}}\\ {{\rm{d}}{e_y}}\\ {{\rm{d}}{e_z}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\rm{d}}{X_0}}\\ {{\rm{d}}{Y_0}}\\ {{\rm{d}}{Z_0}} \end{array}} \right]} \end{array} $$ (2)

      令$\left[ \begin{array}{l} {e_x}\\ {e_y}\\ {e_z} \end{array} \right] = \left[ \begin{array}{l} e_x^i\\ e_y^i\\ e_z^i \end{array} \right] + \left[ \begin{array}{l} {\rm{d}}{e_x}\\ {\rm{d}}{e_y}\\ {\rm{d}}{e_z} \end{array} \right]$,X=[X0  Y0  Z0  μ  α1  α2  α3]T,对式(2)进行整理得:

      $$ {\mathit{\boldsymbol{X}}_2} - {\mathit{\boldsymbol{e}}_2} = \left( {{\mu ^i}{\mathit{\boldsymbol{M}}^i}{\mathit{\boldsymbol{X}}_1} + \Delta \mathit{\boldsymbol{X}}_0^i} \right) + {\mathit{\boldsymbol{A}}^i}{\rm{d}}\mathit{\boldsymbol{X}} - {\mu ^i}{\mathit{\boldsymbol{M}}^i}{\mathit{\boldsymbol{e}}_1} $$ (3)

      式中,

      $$ {\mathit{\boldsymbol{X}}_2} = \left[ {\begin{array}{*{20}{c}} X\\ Y\\ Z \end{array}} \right];{\mathit{\boldsymbol{e}}_2} = \left[ {\begin{array}{*{20}{c}} {{e_X}}\\ {{e_Y}}\\ {{e_Z}} \end{array}} \right];{\mathit{\boldsymbol{X}}_1} = \left[ {\begin{array}{*{20}{c}} x\\ y\\ z \end{array}} \right];\Delta \mathit{\boldsymbol{X}}_0^i = \left[ {\begin{array}{*{20}{c}} {X_0^i}\\ {Y_0^i}\\ {Z_0^i} \end{array}} \right]; $$
      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}^i} = }\\ {\left[ {{\mathit{\boldsymbol{I}}_{3 \times 3}}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{M}}^i}}&{{\mu ^i}\frac{{\partial \mathit{\boldsymbol{M}}}}{{\partial {\alpha _1}}}}&{{\mu ^i}\frac{{\partial \mathit{\boldsymbol{M}}}}{{\partial {\alpha _2}}}}&{{\mu ^i}\frac{{\partial \mathit{\boldsymbol{M}}}}{{\partial {\alpha _3}}}} \end{array}} \right]\left( {{\mathit{\boldsymbol{X}}_1} - \mathit{\boldsymbol{e}}_1^i} \right)} \right];} \end{array} $$
      $$ \begin{array}{l} {\rm{d}}\mathit{\boldsymbol{X}} = {\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{X_0}}&{{\rm{d}}{Y_0}}&{{\rm{d}}{Z_0}}&{{\rm{d}}\mu }&{{\rm{d}}{\alpha _1}}&{{\rm{d}}{\alpha _2}}&{{\rm{d}}{\alpha _3}} \end{array}} \right]^{\rm{T}}};\\ {\mathit{\boldsymbol{e}}_1} = {\left[ {\begin{array}{*{20}{c}} {{e_x}}&{{e_y}}&{{e_z}} \end{array}} \right]^{\rm{T}}}。 \end{array} $$

      采用式(3)构建加权整体最小二乘问题的拉格朗日目标函数Φ

      $$ \begin{array}{l} \mathit{\Phi }\left( {{\mathit{\boldsymbol{e}}_1},{\mathit{\boldsymbol{e}}_2},\mathit{\boldsymbol{X}},\mathit{\boldsymbol{K}}} \right) = \mathit{\boldsymbol{e}}_1^{\rm{T}}\mathit{\boldsymbol{Q}}_1^{ - 1}{\mathit{\boldsymbol{e}}_1} + \mathit{\boldsymbol{e}}_2^{\rm{T}}\mathit{\boldsymbol{Q}}_2^{ - 1}{\mathit{\boldsymbol{e}}_2} + 2{\mathit{\boldsymbol{K}}^{\rm{T}}} \cdot \\ \left( {{\mathit{\boldsymbol{X}}_2} - {\mathit{\boldsymbol{e}}_2} - \left( {{\mu ^i}{\mathit{\boldsymbol{M}}^i}{\mathit{\boldsymbol{X}}_1} + \Delta \mathit{\boldsymbol{X}}_0^i} \right) - {\mathit{\boldsymbol{A}}^i}{\rm{d}}\mathit{\boldsymbol{X}} - {\mu ^i}{\mathit{\boldsymbol{M}}^i}{\mathit{\boldsymbol{e}}_1}} \right) \end{array} $$ (4)

      式中,K为拉格朗日乘子;Q1Q2分别为e1e2的协因数阵。不考虑e1e2的相关性,则观测值的协因数阵可表示为:

      $$ \mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_1}}&{\bf{0}}\\ {\bf{0}}&{{\mathit{\boldsymbol{Q}}_2}} \end{array}} \right] $$ (5)

      通过构建的拉格朗日必要条件推导目标函数的解,对各变量求导并令导数为零:

      $$ \frac{1}{2}\frac{{\partial \mathit{\Phi }}}{{\partial {\mathit{\boldsymbol{e}}_1}}}\left| {_{{{\mathit{\boldsymbol{\tilde e}}}_1},{{\mathit{\boldsymbol{\tilde e}}}_2},{\rm{d}}\mathit{\boldsymbol{\hat X}},\mathit{\boldsymbol{\hat K}}}} \right. = \mathit{\boldsymbol{Q}}_1^{ - 1}{{\mathit{\boldsymbol{\tilde e}}}_1} + {\left( {{\mathit{\boldsymbol{M}}^i}} \right)^{\rm{T}}}{\mu ^i}\mathit{\boldsymbol{\hat K}} = 0 $$ (6)
      $$ \frac{1}{2}\frac{{\partial \mathit{\Phi }}}{{\partial {\mathit{\boldsymbol{e}}_2}}}\left| {_{{{\mathit{\boldsymbol{\tilde e}}}_1},{{\mathit{\boldsymbol{\tilde e}}}_2},{\rm{d}}\mathit{\boldsymbol{\hat X}},\mathit{\boldsymbol{\hat K}}}} \right. = \mathit{\boldsymbol{Q}}_2^{ - 1}{{\mathit{\boldsymbol{\tilde e}}}_2} - \mathit{\boldsymbol{\hat K}} = 0 $$ (7)
      $$ \frac{1}{2}\frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{X}}}}\left| {_{{{\mathit{\boldsymbol{\tilde e}}}_1},{{\mathit{\boldsymbol{\tilde e}}}_2},{\rm{d}}\mathit{\boldsymbol{\hat X}},\mathit{\boldsymbol{\hat K}}}} \right. = - {\left( {{\mathit{\boldsymbol{A}}^i}} \right)^{\rm{T}}}\mathit{\boldsymbol{\hat K}} = 0 $$ (8)
      $$ \begin{array}{*{20}{c}} {\frac{1}{2}\frac{{\partial \mathit{\Phi }}}{{\partial \mathit{\boldsymbol{K}}}}\left| {_{{{\mathit{\boldsymbol{\tilde e}}}_1},{{\mathit{\boldsymbol{\tilde e}}}_2},{\rm{d}}\mathit{\boldsymbol{\hat X}},\mathit{\boldsymbol{\hat K}}}} \right. = {\mathit{\boldsymbol{X}}_2} - {{\mathit{\boldsymbol{\tilde e}}}_2} - \left( {{\mu ^i}{\mathit{\boldsymbol{M}}^i}{\mathit{\boldsymbol{X}}_1} + } \right.}\\ {\left. {\Delta \mathit{\boldsymbol{X}}_0^i} \right) - {\mathit{\boldsymbol{A}}^i}{\rm{d}}\mathit{\boldsymbol{\hat X}} + {\mu ^i}{\mathit{\boldsymbol{M}}^i}{{\mathit{\boldsymbol{\tilde e}}}_1} = 0} \end{array} $$ (9)

      式(6)~(9)中,上标^代表估计值; ~代表预测值。

      由式(6)、(7)和(9)得:

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_2} - {\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{\hat K}} - {\mu ^i}{\mathit{\boldsymbol{M}}^i}{\mathit{\boldsymbol{X}}_1} - \Delta \mathit{\boldsymbol{X}}_0^i - {\mathit{\boldsymbol{A}}^i}{\rm{d}}\mathit{\boldsymbol{\hat X}} - }\\ {{\mu ^i}{\mathit{\boldsymbol{M}}^i}{\mathit{\boldsymbol{Q}}_1}{{\left( {{\mathit{\boldsymbol{M}}^i}} \right)}^{\rm{T}}}{\mu ^i}\mathit{\boldsymbol{\hat K}} = 0} \end{array} $$ (10)

      Li=X2-μiMiX1-ΔX0iQci=Q2+μiMiQ1(Mi)Tμi,由式(10)得:

      $$ \mathit{\boldsymbol{\hat K}} = {\left( {\mathit{\boldsymbol{Q}}_c^i} \right)^{ - 1}}\left( {{\mathit{\boldsymbol{L}}^i} - {\mathit{\boldsymbol{A}}^i}{\rm{d}}\mathit{\boldsymbol{\hat X}}} \right) $$ (11)

      将式(11)代入式(8),得:

      $$ {\left( {{\mathit{\boldsymbol{A}}^i}} \right)^{\rm{T}}}{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)^{ - 1}}\left( {{\mathit{\boldsymbol{L}}^i} - {\mathit{\boldsymbol{A}}^i}{\rm{d}}\mathit{\boldsymbol{\hat X}}} \right) = 0 $$ (12)

      则待求参数的改正数${\rm{d}}\mathit{\boldsymbol{\hat X}}$在第i+1次迭代的估值为:

      $$ {\rm{d}}{{\mathit{\boldsymbol{\hat X}}}^{i + 1}} = {\left( {{{\left( {{\mathit{\boldsymbol{A}}^i}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^i}} \right)^{ - 1}}{\left( {{\mathit{\boldsymbol{A}}^i}} \right)^{\rm{T}}}{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)^{ - 1}}{\mathit{\boldsymbol{L}}^i} $$ (13)

      更新的参数向量为:

      $$ {{\mathit{\boldsymbol{\hat X}}}^{i + 1}} = {{\mathit{\boldsymbol{\hat X}}}^i} + {\rm{d}}{{\mathit{\boldsymbol{\hat X}}}^{i + 1}} $$ (14)

      将式(11)代入式(6)、(7),得到第i+1次迭代的预测残差向量:

      $$ \mathit{\boldsymbol{\tilde e}}_1^{i + 1} = - {\mathit{\boldsymbol{Q}}_1}{\left( {{\mathit{\boldsymbol{M}}^i}} \right)^{\rm{T}}}{\mu ^i}{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)^{ - 1}}\left( {{\mathit{\boldsymbol{L}}^i} - {\mathit{\boldsymbol{A}}^i}{\rm{d}}{{\mathit{\boldsymbol{\hat X}}}^{i + 1}}} \right) $$ (15)
      $$ \mathit{\boldsymbol{\tilde e}}_2^{i + 1} = {\mathit{\boldsymbol{Q}}_2}{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)^{ - 1}}\left( {{\mathit{\boldsymbol{L}}^i} - {\mathit{\boldsymbol{A}}^i}{\rm{d}}{{\mathit{\boldsymbol{\hat X}}}^{i + 1}}} \right) $$ (16)

      将式(14)和式(15)获得的更新后的参数向量和原坐标系统下坐标值的预测残差代入式(13)进行下一次迭代,直到$\left\| {{\rm{d}}{{\mathit{\boldsymbol{\hat X}}}^{i + 1}}} \right\| < {\varepsilon _0}$(ε0为阈值,取ε0=1×10-10),终止迭代,以获得最终的参数估值${\mathit{\boldsymbol{\hat X}}}$和残差项12

      不考虑极小项残差${\mathit{\boldsymbol{A}}^i}{\rm{d}}{{\mathit{\boldsymbol{\hat X}}}^{i + 1}}$,由式(15)和式(16),可得验后单位权方差:

      $$ \sigma _0^2 = \frac{{\mathit{\boldsymbol{\tilde e}}_1^{\rm{T}}\mathit{\boldsymbol{Q}}_1^{ - 1}{{\mathit{\boldsymbol{\tilde e}}}_1} + \mathit{\boldsymbol{\tilde e}}_2^{\rm{T}}\mathit{\boldsymbol{Q}}_2^{ - 1}{{\mathit{\boldsymbol{\tilde e}}}_2}}}{{3n - 7}} = \frac{{{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_c^{ - 1}\mu \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{Q}}_1^{ - 1}{\mathit{\boldsymbol{Q}}_1}{\mathit{\boldsymbol{M}}^{\rm{T}}}\mu \mathit{\boldsymbol{Q}}_c^{ - 1}\mathit{\boldsymbol{L}} + {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_c^{ - 1}{\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{Q}}_2^{ - 1}{\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{Q}}_c^{ - 1}\mathit{\boldsymbol{L}}}}{{3n - 7}} = \frac{{{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_c^{ - 1}\mathit{\boldsymbol{L}}}}{{3n - 7}} $$ (17)

      式中,n为参与计算的公共点位个数。由于迭代过程完成,所有符号的右上标去除。

    • 基于等价权抗差原理进行参数的抗差估计,需要构建检验量,其关键在于残差向量的协因数阵的推导。

      由式(15)和式(16)得到第i+1次迭代的预测残差向量:

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}^{i + 1}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tilde e}}_1^{i + 1}}\\ {\mathit{\boldsymbol{\tilde e}}_2^{i + 1}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - {\mathit{\boldsymbol{Q}}_1}{{\left( {{\mathit{\boldsymbol{M}}^i}} \right)}^{\rm{T}}}{\mu ^i}{{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)}^{ - 1}}}\\ {{\mathit{\boldsymbol{Q}}_2}{{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)}^{ - 1}}} \end{array}} \right] \cdot }\\ {\left( {{\mathit{\boldsymbol{L}}^i} - {\mathit{\boldsymbol{A}}^i}{\rm{d}}{{\mathit{\boldsymbol{\hat X}}}^{i + 1}}} \right) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{C}}^i}}\\ {{\mathit{\boldsymbol{D}}^i}} \end{array}} \right] \cdot {\mathit{\boldsymbol{R}}^{i + 1}}} \end{array} $$ (18)

      式中,为了方便起见,用CDR代替对应位置的矩阵。则第i+1次迭代的预测残差向量的协因数阵为:

      $$ \mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{V}}^{i + 1} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{C}}^i}}\\ {{\mathit{\boldsymbol{D}}^i}} \end{array}} \right] \cdot {\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{R}}^{i + 1}}}} \cdot \left[ {\begin{array}{*{20}{c}} {{{\left( {{\mathit{\boldsymbol{C}}^i}} \right)}^{\rm{T}}}}&{{{\left( {{\mathit{\boldsymbol{D}}^i}} \right)}^{\rm{T}}}} \end{array}} \right] $$ (19)

      由式(13)和(18)可得:

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}^{i + 1}} = {\mathit{\boldsymbol{L}}^i} - {\mathit{\boldsymbol{A}}^i}{\rm{d}}{{\mathit{\boldsymbol{\hat X}}}^{i + 1}} = {\mathit{\boldsymbol{L}}^i} - {\mathit{\boldsymbol{A}}^i} \cdot {{\left( {{{\left( {{\mathit{\boldsymbol{A}}^i}} \right)}^{\rm{T}}} \cdot {{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)}^{ - 1}} \cdot {\mathit{\boldsymbol{A}}^i}} \right)}^{ - 1}} \cdot {{\left( {{\mathit{\boldsymbol{A}}^i}} \right)}^{\rm{T}}} \cdot {{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)}^{ - 1}} \cdot {\mathit{\boldsymbol{L}}^i} = }\\ {\left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{A}}^i} \cdot {{\left( {{{\left( {{\mathit{\boldsymbol{A}}^i}} \right)}^{\rm{T}}} \cdot {{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)}^{ - 1}} \cdot {\mathit{\boldsymbol{A}}^i}} \right)}^{ - 1}} \cdot {{\left( {{\mathit{\boldsymbol{A}}^i}} \right)}^{\rm{T}}} \cdot {{\left( {\mathit{\boldsymbol{Q}}_c^i} \right)}^{ - 1}}} \right) \cdot {\mathit{\boldsymbol{L}}^i} = \left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{F}}^i}} \right) \cdot {\mathit{\boldsymbol{L}}^i}} \end{array} $$ (20)

      式中,用F代替对应位置的矩阵。则Ri+1的协因数阵为:

      $$ {\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{R}}^{i + 1}}}} = \left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{F}}^i}} \right) \cdot {\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{L}}^i}}} \cdot {\left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{F}}^i}} \right)^{\rm{T}}} $$ (21)

      由式(10)有:

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}^i} = {\mathit{\boldsymbol{X}}_2} - {\mu ^i}{\mathit{\boldsymbol{M}}^i}{\mathit{\boldsymbol{X}}_1} - \Delta \mathit{\boldsymbol{X}}_0^i = }\\ {\left[ {\begin{array}{*{20}{c}} { - {\mu ^i}{\mathit{\boldsymbol{M}}^i}}&\mathit{\boldsymbol{I}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_1}}\\ {{\mathit{\boldsymbol{X}}_2}} \end{array}} \right] - \Delta \mathit{\boldsymbol{X}}_0^i} \end{array} $$ (22)

      Li的协因数阵为:

      $$ {\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{L}}^i}}} = \left[ {\begin{array}{*{20}{c}} { - {\mu ^i}{\mathit{\boldsymbol{M}}^i}}&\mathit{\boldsymbol{I}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_1}}&{\bf{0}}\\ {\bf{0}}&{{\mathit{\boldsymbol{Q}}_2}} \end{array}} \right]{\left[ {\begin{array}{*{20}{c}} { - {\mu ^i}{\mathit{\boldsymbol{M}}^i}}&\mathit{\boldsymbol{I}} \end{array}} \right]^{\rm{T}}} $$ (23)

      将式(23)代入式(21),并将式(21)代入式(19),即得到预测残差协因数阵的表达式。由于上述推导过程获得的预测残差协因数阵与迭代过程中的参数和系数阵残差有关,故预测残差协因数阵需要在迭代中不断更新。

      在此基础上,对式(18)中的预测残差进行标准化:

      $$ {{\tilde v}_j} = \frac{{{v_j}}}{{{\sigma _0}\left( {\sqrt {{q_{{v_j}}}} } \right)}}\;\;\;\left( {{q_{{v_j}}} \ne 0} \right) $$ (24)

      式中,j为标准化残差向量中的第j个元素;vj表示预测残差向量V中的第j个元素;σ0为单位权中误差;qvj为残差协因数阵QV中对角线上第j个元素。σ0由中位函数计算,可表示为:

      $$ {\sigma _0} = \mathop {{\rm{med}}}\limits_{j = 1}^n \left( {\left| {{v_j}/\sqrt {{q_{{v_j}}}} } \right|} \right) \cdot 1.4826\left( {{q_{{v_j}}} \ne 0} \right) $$ (25)

      因为标准化残差进行抗差估计可选的权因子函数众多[27],本文采用IGG Ⅲ权因子函数,为了计算方便,采用其对应协因数因子函数Rjj[26]

      $$ {R_{jj}} = \left\{ \begin{array}{l} 1.0,\left| {{{\tilde v}_j}} \right| \le {k_0}\\ \frac{{\left| {{{\tilde v}_j}} \right|}}{{{k_0}}}{\left( {\frac{{{k_1} - {k_0}}}{{{k_1} - \left| {{{\tilde v}_j}} \right|}}} \right)^2},{k_0} < \left| {{{\tilde v}_j}} \right| \le {k_1}\\ 1 \times {10^{10}},{k_1} < \left| {{{\tilde v}_j}} \right| \end{array} \right. $$ (26)

      式中,第3段理论上应该为无穷大,为了方便计算,可采用极大值替代(如1×1010),可满足实际要求;k0取为2.0~3.0;k1取为4.5~8.5[27]。采用双因子模型,由式(26)可得观测向量和系数矩阵中观测值的等价协因数:

      $$ {{\bar q}_{{j_k}}} = {q_{{j_k}}}\sqrt {{R_{jj}}} \sqrt {{R_{kk}}} $$ (27)

      式中,qjkqjk分别为观测值的等价协因数阵Q和协因数阵Q中对应位置元素。采用Q替代Q,构建类似于式(4)的加权整体最小二乘问题的拉格朗日目标函数,对各变量求导并令导数为零,类似于式(6)~(9),可得第i+1次迭代的待求参数抗差估计的改正数为:

      $$ {\rm{d}}{{\mathit{\boldsymbol{\hat {\bar X}}}}^{i + 1}} = {\left( {{{\left( {{\mathit{\boldsymbol{A}}^i}} \right)}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{\bar Q}}_c^i} \right)}^{ - 1}}{\mathit{\boldsymbol{A}}^i}} \right)^{ - 1}}{\left( {{\mathit{\boldsymbol{A}}^i}} \right)^{\rm{T}}}{\left( {\mathit{\boldsymbol{\bar Q}}_c^i} \right)^{ - 1}}{\mathit{\boldsymbol{L}}^i} $$ (28)

      式中,Qci=Q2+μiMiQ1(Mi)Tμi。采用类似于式(13)~(15)的迭代过程,在给定的阈值条件下,可获得参数的最终估计值及其精度评定指标。

    • 设有均匀分布的25个点,原坐标系下坐标为(x1y1z1);坐标转换的7参数分别为平移参数(X0Y0Z0)=(1 000,1 000,1 000)m;尺度因子μ=2.0;旋转参数(α1α2α3)=(1.0,0.5,1.5)rad;计算得到目标坐标系下的坐标为(x2y2z2)。

      随机选择18个点作为公共点,另外7个点作为检验点。采用正态分布函数,对公共点在原坐标系和目标坐标系下的坐标真值加入不等精度的随机误差,其数学期望为0,标准差在0~0.05 m间随机产生,表示为(dx1,dy1,dz1)、(dx2,dy2,dz2), 进而获得18个点在两套坐标系下的仿真坐标(x1+dx1y1+dy1z1+dz1)、(x2+dx2y2+dy2z2+dz2)。粗差的位置随机产生,大小为标准差的5~20倍,加入粗差个数分别为1、3、5。同样,对7个检验点在原坐标系下的坐标真值加入随机误差,随机误差的产生方式与公共点相同,并保留其在目标坐标系下的坐标真值,用于分析不同方案的性能;仿真次数为500。

      采用6种方案进行参数估计, 无粗差情况下:(1)方案1, 经典最小二乘算法(LS);(2)方案2, 文献[14]提出的通用整体最小二乘算法(WTLS);(3)方案3, 本文提出的基于GH模型的整体最小二乘算法(GH-WTLS)。有粗差情况下:(4)方案4,本文提出的基于GH模型的整体最小二乘算法(GH-WTLS);(5)方案5, 采用文献[24]提出的RWTLS方法对GH模型进行参数估计(RWTLS);(6)方案6, 本文提出的GH模型的抗差估计算法(GH-RWTLS)。

      表 1~3分别给出了加入不同数目粗差、6种方案求得的参数估值的均方根误差。可知:(1)方案2和3的均方根误差明显小于方案1。即不含有粗差时,WTLS方法的估计结果优于LS,说明了整体最小二乘顾及系数矩阵中随机误差的合理性;(2)方案2与方案3的各项指标一致,说明了本文方法的正确性,且本文方法较文献[14]的方法,估计参数数目由61个减少至7个,矩阵计算维数大大降低;(3)比较方案3与方案4,受粗差影响,方案4参数估值的均方根误差明显增大,即整体最小二乘算法抵抗粗差的能力较差;(4)方案5、6与方案4相比较,参数估值的均方根误差得到了不同程度的减小,即采用抗差估计方法可以削弱粗差对参数估值的影响,提高估计精度;(5)无论加入几个粗差,方案6估值的精确度均优于方案5,且随着粗差的增多,优势更为明显。主要是由于在非等精度观测值条件下,方案5构建的检验量未顾及结构空间的影响。进一步比较方案5、6,图 1给出了加入3个粗差时,方案5、6的参数估计值与真值的差值序列,方案6对应的序列收敛程度明显优于方案5,参数估值更为稳健。

      表 1  存在1个粗差时各方案的均方根误差对比

      Table 1.  Root Mean Square Errors of Estimated Parameters for Different Schemes with One Outlier

      方案X0 /mY0 /mZ0 /mμ/(10-5 rad)α1 /(10-5 rad)α2 /(10-5 rad)α3 /(10-5 rad)
      方案119.302 110.579 615.580 12.6001.5101.1601.750
      方案23.922 12.666 33.613 20.6320.4050.3510.366
      方案33.92 212.666 33.613 20.6320.4050.3510.366
      方案46.470 33.981 95.938 81.0100.6440.5570.570
      方案54.502 93.163 64.381 60.6960.4530.4100.452
      方案64.194 92.907 24.106 60.6660.4350.3820.407

      表 2  存在3个粗差时各方案的均方根误差对比

      Table 2.  Root Mean Square Errors of Estimated Parameters for Different Schemes with Three Outliers

      方案X0 /mY0 /mZ0 /mμ/(10-5 rad)α1 /(10-5 rad)α2 /(10-5 rad)α3 /(10-5 rad)
      方案110.307 912.340 718.474 01.4601.3901.7301.130
      方案23.752 92.956 42.844 90.5070.2890.3110.417
      方案33.752 92.956 42.844 90.5070.2890.3110.417
      方案48.797 86.880 56.709 91.2300.6640.7700.952
      方案57.068 05.744 25.810 60.9780.5920.6280.798
      方案64.695 83.896 13.931 20.6710.4040.4130.518

      表 3  存在5个粗差时各方案的均方根误差对比

      Table 3.  Root Mean Square Errors of Estimated Parameters for Different Schemes with Five Outliers

      方案X0 /mY0 /mZ0 /mμ/(10-5 rad)α1 /(10-5 rad)α2 /(10-5 rad)α3 /(10-5 rad)
      方案112.748 611.219 615.764 92.2501.7901.1501.400
      方案23.414 03.380 14.576 30.6100.4690.3800.424
      方案33.414 03.380 14.576 30.6100.4690.3800.424
      方案410.173 610.600 614.184 51.9601.4801.1201.320
      方案58.900 88.467 611.919 81.5801.2900.9561.110
      方案66.064 66.039 37.999 31.0900.8180.6390.757

      图  1  方案5和6估值与真值的差值序列(3个粗差)

      Figure 1.  Difference Values Between Truth Values and Estimated Results of Schemes 5 and 6 with Three Outliers

      为了检验上述不同粗差个数条件下、不同方案求得的转换参数的外符合精度,对7个检查点(无粗差)进行坐标转换,并与对应的目标坐标下的坐标真值进行对比。重点比较方案5和方案6的转换结果,表 4~6分别给出了采用在1个粗差、3个粗差和5个粗差情况下求得的转换参数,进行检查点坐标转换,以及在XYZ 3个方向上的均方根误差。总体上,对于不同的粗差数目,采用方案6得到的转换参数进行坐标转换,结果的均方根误差小于方案5;且随着粗差数目的增加,方案6的优势越显著,这与表 1~3结果基本一致,表明了方案6的内外符合性均优于方案5,具有更好的抗差性能。

      表 4  方案5和方案6计算检查点的均方根误差(1个粗差)/m

      Table 4.  Root Mean Square Errors of Check Points for Schemes 5 and 6 with One Outlier/m

      方向方案点位1点位2点位3点位4点位5点位6点位7
      X方案50.079 20.081 30.072 90.077 60.081 50.080 30.080 3
      方案60.078 40.080 30.072 40.077 10.080 80.079 00.078 7
      Y方案50.096 20.040 70.076 90.039 80.065 60.068 20.041 0
      方案60.095 90.039 10.077 00.038 60.065 30.068 20.040 1
      Z方案50.084 80.052 20.067 50.059 80.077 60.082 80.057 3
      方案60.084 70.052 20.066 80.058 30.076 40.081 90.056 1

      表 5  方案5和方案6计算检查点的均方根误差(3个粗差)/m

      Table 5.  Root Mean Square Errors of Check Points for Schemes 5 and 6 with Three Outliers/m

      方向方案点位1点位2点位3点位4点位5点位6点位7
      X方案50.078 50.081 10.073 90.082 20.082 70.080 30.082 9
      方案60.077 80.078 50.071 70.078 20.080 40.078 50.078 0
      Y方案50.098 70.048 10.082 70.044 60.065 10.073 00.044 9
      方案60.094 60.041 90.078 40.038 70.063 90.070 70.039 2
      Z方案50.089 30.056 90.068 30.061 80.076 80.083 40.059 3
      方案60.087 10.052 90.066 20.057 60.076 20.082 10.057 2

      表 6  方案5和方案6计算检查点的均方根误差(5个粗差)/m

      Table 6.  Root Mean Square Errors of Check Points for Schemes 5 and 6 with Five Outliers/m

      方向方案点位1点位2点位3点位4点位5点位6点位7
      X方案50.084 40.09 220.080 60.089 30.087 30.090 70.097 6
      方案60.079 70.083 50.075 00.082 10.081 30.081 60.085 5
      Y方案50.104 00.057 70.087 70.052 20.072 80.079 40.057 6
      方案60.097 90.047 50.083 10.041 90.069 50.072 10.046 8
      Z方案50.099 50.066 10.078 10.073 40.085 40.087 90.068 1
      方案60.091 90.060 10.070 80.062 10.078 90.083 50.060 0
    • 对三维坐标转换的GH模型,本文采用Newton-Gauss迭代算法,推导了该模型解算方法。在此基础上,基于标准化残差构造权因子函数,给出了该模型的抗差解法。通过理论推导和仿真实验分析,可得出:(1)GH模型不受限制条件影响,适用于大角度的三维坐标转换问题,且相较于文献[14],解算结果完全一致,而待估参数的数目大大降低,其计算效率有一定程度的提高;(2)本文提出的抗差解法综合考虑了非等精度条件下可能出现的粗差对观测空间与结构空间的影响,仿真结果表现出了良好的稳健性。

参考文献 (27)

目录

    /

    返回文章
    返回