﻿ 三维坐标转换的稳健加权总体最小二乘算法
 文章快速检索 高级检索
 大地测量与地球动力学  2020, Vol. 40 Issue (10): 1027-1033  DOI: 10.14075/j.jgg.2020.10.007

### 引用本文

WANG Leyang, XU Ranran. Robust Weighted Total Least Squares Algorithm for Three-Dimensional Coordinate Transformation[J]. Journal of Geodesy and Geodynamics, 2020, 40(10): 1027-1033.

### Foundation support

National Natural Science Foundation of China, No. 41874001, 41664001;Innovation Fund Designated for Graduate Students of ECUT, No.DHYC-201925.

### 第一作者简介

WANG Leyang, PhD, associate professor, majors in geodetic inversion and geodetic data processing, E-mail:wleyang@163.com.

### 文章历史

1. 东华理工大学测绘工程学院，南昌市广兰大道418号，330013

1 基于四元数的加权对称相似变换(WSSTUQ)

 $\left[ {\begin{array}{*{20}{l}} {X + {e_X}}\\ {Y + {e_Y}}\\ {Z + {e_Z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{T_X}}\\ {{T_Y}}\\ {{T_Z}} \end{array}} \right] + k\mathit{\boldsymbol{R}}\left[ {\begin{array}{*{20}{l}} {x + {e_x}}\\ {y + {e_y}}\\ {z + {e_z}} \end{array}} \right]$ (1)

 $\psi ({\mathit{\boldsymbol{e}}_{{1}}},{\mathit{\boldsymbol{e}}_{{2}}},\mathit{\boldsymbol{\xi }}) = h\left( {{\mathit{\boldsymbol{y}}_{{1}}} + {\mathit{\boldsymbol{e}}_{{1}}},{\mathit{\boldsymbol{y}}_{{2}}} + {\mathit{\boldsymbol{e}}_{{2}}},\mathit{\boldsymbol{\xi }}} \right) = 0$ (2)

 $\begin{array}{l} f({\mathit{\boldsymbol{e}}_{{1}}},{\mathit{\boldsymbol{e}}_{{2}}},\mathit{\boldsymbol{\xi }}) \approx {\mathit{\boldsymbol{B}}_{{{{1}}_{{0}}}}}({\mathit{\boldsymbol{e}}_{{1}}} - \mathit{\boldsymbol{e}}_{{1}}^{{0}}) + {\mathit{\boldsymbol{B}}_{{{{2}}_{{0}}}}}({\mathit{\boldsymbol{e}}_{{2}}} - \mathit{\boldsymbol{e}}_{{2}}^{{0}}) + \\ \mathit{\boldsymbol{A}}(\mathit{\boldsymbol{\xi }} - {\mathit{\boldsymbol{\xi }}^{{0}}}) + \psi (\mathit{\boldsymbol{e}}_{{1}}^{{0}},\mathit{\boldsymbol{e}}_{{2}}^{{0}},{\mathit{\boldsymbol{\xi }}^{{0}}}) = {{0}} \end{array}$ (3)

 $\mathit{\Phi} = \mathit{\boldsymbol{e}}_{{1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{1}}^{ - {\rm{1}}}{\mathit{\boldsymbol{e}}_{{1}}} + \mathit{\boldsymbol{e}}_{{2}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{2}}^{ - {\rm{1}}}{\mathit{\boldsymbol{e}}_{{2}}} = \min$ (4)

 $\begin{array}{l} \mathit{\Phi} ({\mathit{\boldsymbol{e}}_{{1}}},{\mathit{\boldsymbol{e}}_{{2}}},\mathit{\boldsymbol{\xi }},\mathit{\boldsymbol{\lambda }}) = \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{\lambda }}^{\rm{T}}}\psi ({\mathit{\boldsymbol{e}}_{{1}}},{\mathit{\boldsymbol{e}}_{{2}}},\mathit{\boldsymbol{\xi }}) = \min \end{array}$ (5)

 ${\left. {\frac{{\partial \mathit{\Phi} }}{{2\partial {\mathit{\boldsymbol{e}}_{{1}}}}}} \right|_{{{{e}}_{{1}}}{{,}}{{{e}}_{{2}}}{{, {{ ξ}} , {{ λ}} }}}} = \mathit{\boldsymbol{Q}}_{{1}}^{ - {\rm{1}}}{\mathit{\boldsymbol{e}}_{{1}}} - \mathit{\boldsymbol{B}}_{{{{1}}_{{0}}}}^{\rm{T}}\mathit{\boldsymbol{\lambda }} = {{0}}$ (6)
 ${\left. {\frac{{\partial \mathit{\Phi} }}{{2\partial {\mathit{\boldsymbol{e}}_{{2}}}}}} \right|_{{{{e}}_{{1}}}{{,}}{{{e}}_{{2}}}{{, {{ ξ}} , {{ λ}} }}}} = \mathit{\boldsymbol{Q}}_{{2}}^{ - {\rm{1}}}{\mathit{\boldsymbol{e}}_{{2}}} - \mathit{\boldsymbol{B}}_{{{{2}}_{{0}}}}^{\rm{T}}\mathit{\boldsymbol{\lambda }} = {{0}}$ (7)
 ${\left. {\frac{{\partial \mathit{\Phi} }}{{2\partial \mathit{\boldsymbol{\xi }}}}} \right|_{{{{e}}_{{1}}}{{,}}{{{e}}_{{2}}}{{, {{ ξ}} , {{ λ}} }}}} = {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{\lambda }} = {{0}}$ (8)
 ${\left. {\frac{{\partial \mathit{\Phi} }}{{2\partial \mathit{\boldsymbol{\lambda }}}}} \right|_{{{{e}}_{{1}}}{{,}}{{{e}}_{{2}}}{{, {{ ξ}} , {{ λ}} }}}} = \psi ({\mathit{\boldsymbol{e}}_{{1}}},{\mathit{\boldsymbol{e}}_{{2}}},\mathit{\boldsymbol{\zeta }}) = {{0}}$ (9)

 $\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{1}}} + {\mathit{\boldsymbol{B}}_{{{{2}}_{{0}}}}}{\mathit{\boldsymbol{Q}}_{{2}}}\mathit{\boldsymbol{B}}_{{{{2}}_{{0}}}}^{\rm{T}}}&\mathit{\boldsymbol{A}}\\ {{\mathit{\boldsymbol{A}}^{\rm{T}}}}&{{0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat \lambda }}}^{{1}}}}\\ {{{\mathit{\boldsymbol{\hat \xi }}}^{{1}}} - {\mathit{\boldsymbol{\xi }}^{{0}}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{w}}\\ {{0}} \end{array}} \right] = {{0}}$ (10)

 $\left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat e}}}_{{1}}}}\\ {{{\mathit{\boldsymbol{\hat e}}}_{{2}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{1}}}}\\ {{\mathit{\boldsymbol{Q}}_{{2}}}\mathit{\boldsymbol{B}}_{{{{2}}_{{0}}}}^{{\rm T}}} \end{array}} \right]{\mathit{\boldsymbol{\hat \lambda }}^{{1}}}$ (11)

 $\sigma _0^2 = \frac{{\mathit{\boldsymbol{e}}_{{1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{1}}^{ - {\rm{1}}}{\mathit{\boldsymbol{e}}_{{1}}} + \mathit{\boldsymbol{e}}_{{2}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{2}}^{ - {\rm{1}}}{\mathit{\boldsymbol{e}}_{{2}}}}}{{3m - 7}}$ (12)

 ${\rm{bias}}{}_{{ {{ ξ}} }} = \sum\nolimits\limits_{i = 1}^n {({\mathit{\boldsymbol{\xi }}_{{i}}} - {\mathit{\boldsymbol{\xi }}_{{\rm{tv}}}})} /n$ (13)
 ${\rm{RMSE}} = \sqrt {\sum\nolimits\limits_{i = 1}^n {{{({\mathit{\boldsymbol{\xi }}_{{i}}} - {\mathit{\boldsymbol{\xi }}_{{\rm{tv}}}})}^2}/n} }$ (14)

2 基于四元数的加权对称相似变换稳健估计 2.1 稳健加权对称相似变换(RWSST)

 ${\mathit{\boldsymbol{V}}_{{e}}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat e}}_{{1}}^{{\rm{(}}i + {\rm{1)}}}}\\ {\mathit{\boldsymbol{\hat e}}_{{2}}^{{\rm{(}}i + {\rm{1)}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{M}}\\ \mathit{\boldsymbol{N}} \end{array}} \right]\mathit{\boldsymbol{R}}$ (15)

 ${\mathit{\boldsymbol{Q}}_{{{{V}}_{{e}}}}} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{Q}}_{{R}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}}&{\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{Q}}_{{R}}}{\mathit{\boldsymbol{N}}^{\rm{T}}}}\\ {\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{Q}}_{{R}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}}&{\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{Q}}_{{R}}}{\mathit{\boldsymbol{N}}^{\rm{T}}}} \end{array}} \right]$ (16)

R=LA(i) 1(i+1)=LA(i)((AT)(i)(Qc(i))－1 (A(i)))－1(AT)(i)(Qc(i))－1L= (I-A(i)((AT)(i) (Qc(i))－1(A(i)))－1(AT)(i)(Qc(i))－1)L=(IF(i))L式中，F=A(i)((AT)(i)(Qc(i))－1(A(i)))－1(AT)(i) (Qc(i))－1

R的协因数阵为：

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

 $\begin{array}{l} \mathit{\boldsymbol{L}} = {\mathit{\boldsymbol{y}}_{{1}}} - {\left[ {\begin{array}{*{20}{c}} {{T_X}}&{{T_Y}}&{{T_Z}} \end{array}} \right]^{\rm{T}}} - k\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{y}}_{{2}}} = \\ \;\;\;\;\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{{3}}}}&{ - k\mathit{\boldsymbol{R}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{y}}_{{1}}}}\\ {{\mathit{\boldsymbol{y}}_{{2}}}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{T_X}}\\ {{T_Y}}\\ {{T_Z}} \end{array}} \right] \end{array}$ (18)

 ${\mathit{\boldsymbol{Q}}_{{L}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{{3}}}}&{ - k\mathit{\boldsymbol{R}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{{1}}}}&{{0}}\\ {{0}}&{{\mathit{\boldsymbol{Q}}_{{2}}}} \end{array}} \right]{\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{{3}}}}&{ - k\mathit{\boldsymbol{R}}} \end{array}} \right]^{\rm{T}}}$ (19)

 ${{{\bar v}}_{{{eg}}}} = \frac{{{{{v}}_{{{eg}}}}}}{{{\sigma _0}\sqrt {{{{Q}}_{{{{v}}_{{{eg}}}}}}} }}$ (20)

 ${\sigma _0} = 1.483\;\mathop {\rm med}\limits_{g = 1}^n \left( {\left| {\frac{{{v_{eg}}}}{{\sqrt {{Q_{{v_{eg}}}}} }}} \right|} \right)$ (21)

 ${r_{gg}} = \left\{ {\begin{array}{*{20}{l}} {1.0,\left| {{{\bar v}_{eg}}} \right| < {k_0}}\\ {\exp \left( {1 - {{\left( {\frac{{{{\bar v}_{eg}}}}{{{k_1}}}} \right)}^2}} \right),{k_0} \le \left| {{{\bar v}_{eg}}} \right| < {k_1}}\\ {{{10}^{ - 10}},\left| {{{\bar v}_{eg}}} \right| \ge {k_1}} \end{array}} \right.$ (22)

 ${\bar p_g} = {p_g}{r_{gg}}$ (23)

 ${\bar p_{gl}} = {p_{gl}}\sqrt {{r_{gg}}{r_{ll}}}$ (24)

2.2 稳健加权对称相似变换(RWSST)迭代过程

1) 给定七参数初值ξ0，一般取k=1，其余为0；给定观测向量误差初值e10e20，一般取误差初值为0，输入观测值权阵P

2) 将输入初值代入式(6)~式(9)，计算矩阵B1B2A

3) 通过式(10)求解转换参数改正数dξ(i+1)，通过式(11)求解预测残差Ve(i+1)

4) 将QR=(IF(i))QL(IF(i))T代入式(16)，求解预测残差协因数阵Q(i+1)Ve

5) 如果$\left\| {{\rm{d}}{{\mathit{\boldsymbol{\xi }}}^{(i + 1)}}} \right\| = \left\| {{{\widehat {\mathit{\boldsymbol{\xi }}}}^{(i + 1)}} - {{\widehat {\mathit{\boldsymbol{\xi }}}}^{(i)}}} \right\| > \varepsilon$，则${{\mathit{\boldsymbol{\hat \xi }}}^{(i)}} = {{\mathit{\boldsymbol{\hat \xi }}}^{(i + 1)}},{\mathit{\boldsymbol{V}}}_e^{(i)} = {\mathit{\boldsymbol{V}}}_e^{(i + 1)}$，重复步骤2)~4)；否则，迭代终止；

6) 将获得的残差Ve及残差协因数阵QVe代入式(21)，求解抗差性的单位权中误差σ0，然后将σ0代入式(20)构造标准化残差veg，最后通过veg构造权因子函数rgg

7) 通过式(23)或式(24)计算pgpgl，组建等价权矩阵p，得到观测值的新权值；

8) 重复步骤1)~5)，最终得到七参数估值及观测值残差；

9) 由式(12)进行精度评定，输出σ02

3 模拟算例

 图 1 坐标分布 Fig. 1 Distribution of coordinates

 图 2 含有6个粗差点条件下采用方案3和方案4所得参数偏差二范数对比 Fig. 2 The norm comparison of two parameter deviation when schemes 3 and 4 contains six gross errors

1) 当观测值只含有随机误差时方案1的结果最优；当观测值同时含有少量粗差时，方案2所得参数估值偏离真值；随着所含粗差点个数的增加，方案2所求参数估值偏离真值的速度加快，其结果严重偏离真值。这是因为方案2只能处理观测值含有随机误差的情况，无法顾及含有粗差的情况，不具备抵抗粗差的能力。

2) 表 1~3显示，当粗差点个数在1~3范围内时，方案3的IGG方法与方案4的稳健加权对称相似变换算法求解的转换参数均方根误差与方案1均方根误差相差甚微，即3个平移参数均方根误差差值均在5 mm以内，尺度比例因子均方根误差差值在10-7量级左右，由于3个欧拉角参数均方根误差差值变动较小，此处不作比较。因此，在含有较少粗差的情况下，方案4优于方案3，虽不明显但仍能证明方案4算法更合理，所求参数的解更精确可靠。

3) 表 4~6显示，当含有4~6个粗差点时，方案4所得七参数的均方根误差明显低于方案3，且在含有6个粗差点时，方案3七参数偏差的二范数明显比方案4的大，方案4七参数偏差的二范数分布更集中。随着粗差点个数的增多，方案3与方案4的七参数均方根误差增大速度加快，二者抵抗粗差的能力明显减弱，但方案4所得的各参数均方根误差仍小于方案3，且方案4所得的七参数均方根误差变化不明显，因此方案4明显优于方案3。

 图 3 当目标坐标系下第1、第16、第20号点含有粗差时采用稳健加权对称相似变换迭代后的权因子 Fig. 3 The weight factor after iteration of the RWSST method while 1 st, 16 th and 20 th points in target system include the gross errors

 图 4 当目标坐标系下第1、第16、第20、第35、第41号点含有粗差时采用IGG方法迭代后的权因子 Fig. 4 The weight factor after iteration of IGG method while 1 st, 16 th, 20 th, 35 th and 41 th points in target system include the gross errors

 图 5 当目标坐标系下第1、第16、第20、第35、第41号点含有粗差时采用稳健加权总体最小二乘方法迭代后的权因子 Fig. 5 The weight factor after iteration of the RWSST method while 1 st, 16 th, 20 th, 35 th and 41 th points in target system include the gross errors

4 实例

5 结语

 [1] 王乐洋, 余航. 总体最小二乘方法的适用性研究[J]. 大地测量与地球动力学, 2014, 34(3): 121-124 (Wang Leyang, Yu Hang. Study on Applicability of Total Least Squares Method in Surveying Adjustment[J]. Journal of Geodesy and Geodynamics, 2014, 34(3): 121-124) (0) [2] 王乐洋, 于冬冬. 病态总体最小二乘问题的虚拟观测解法[J]. 测绘学报, 2014(6): 575-581 (Wang Leyang, Yu Dongdong. Virtual Observation Method to Ill-Posed Total Least Squares Problem[J]. Acta Geodaetica et Cartographica Sinica, 2014(6): 575-581) (0) [3] Xu P L, Liu J N. 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 (0) [4] 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 (0) [5] Mercan H, Akyilmaz O, Aydin C. Solution of the Weighted Symmetric Similarity Transformations Based on Quaternions[J]. Journal of Geodesy, 2018, 92(10): 1113-1130 DOI:10.1007/s00190-017-1104-0 (0) [6] 王新洲. 高等测量平差[M]. 北京: 测绘出版社, 2006 (Wang Xinzhou. Higher Surveying Adjustment[M]. Beijing: Surveying and Mapping Press, 2006) (0) [7] 周江文, 黄幼才, 杨元喜, 等. 抗差最小二乘法[M]. 武汉: 华中理工大学出版社, 1997 (Zhou Jiangwen, Huang Youcai, Yang Yuanxi, et al. Robustified Least Squares Approaches[M]. Wuhan: Huazhong University of Science and Engineering Press, 1997) (0) [8] 王乐洋, 鲁铁定. 总体最小二乘平差法的误差传播定律[J]. 大地测量与地球动力学, 2014, 34(2): 55-59 (Wang Leyang, Lu Tieding. Propagation Law of Errors in Total Least Squares Adjustment[J]. Journal of Geodesy and Geodynamics, 2014, 34(2): 55-59) (0) [9] Yang Y, Song L, Xu T H. Robust Estimator for Correlated Observations Based on Bifactor Equivalent Weights[J]. Journal of Geodesy, 2002, 76(6): 353-358 (0) [10] 刘超, 王彬, 赵兴旺, 等. 三维坐标转换的高斯-赫尔默特模型及其抗差解法[J]. 武汉大学学报:信息科学版, 2018, 43(9): 1320-1327 (Liu Chao, Wang Bin, Zhao Xingwang, et al. 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) (0)
Robust Weighted Total Least Squares Algorithm for Three-Dimensional Coordinate Transformation
WANG Leyang1     XU Ranran1
1. Faculty of Geomatics, East China University of Technology, 418 Guanglan Road, Nanchang 330013, China
Abstract: Based on the existing weighted symmetric similarity transformation, only the case where the observation value contains random error is considered, and the case where the observation value contains the gross error is not considered. This paper further verifies that the weighted symmetric similarity transformation is not robust. Based on the weighted symmetric similarity transformation, the method of selecting weight iteration is adopted to make the robust weighted symmetric similarity transformation. The model obtains unit weight mean square error with robustness by the median method and utilizes the standardized residual to construct the weight factor function, to obtain a reliable parametric solution. Comparative analysis shows that: when the observations contain 4-6 gross errors, the method of this paper can be used to detect more data that may have gross errors, the weighting factors given are more reasonable, and the obtained parameter solution has higher accuracy and stronger stability.
Key words: coordinate transformation; total least squares; robust estimation; median method