加权整体最小二乘算法的改进*

杨仕平 范东明 龙玉春

1)西南交通大学地球科学与环境工程学院测绘工程系，成都 610031

2)重庆市合川龙市中学，合川401564

1 引言

从Golub 和van Loan［1］提出利用基于奇异值分解法(SVD，Singular Value Decom-position)的整体最小二乘法求解EIV 模型至今，整体最小二乘法一直在不断地发展。鉴于SVD 方法计算复杂不利于编程，研究人员陆续提出了整体最小二乘的其他解法，如迭代算法、完全正交分解法等改进方法［2－4］。

在大地测量和工程测量数据处理中，绝大部分都可以归结为参数估计问题，这些参数估计问题大致可以分为普通最小二乘估计问题和整体最小二乘估计问题，而整体最小二乘估计问题可进一步分为系数阵的元素全部含有随机误差及只有部分含有随机误差两种情况。

针对测绘数据处理中系数矩阵有特殊结构的EIV 模型，陆钰等［5］提出用LS-TLS 将系数矩阵中不需要修正的列固定，但它不能固定所有的不需要修正的元素，不能考虑重复元素;Schaffrin 等［6］提出在加权整体最小二乘法解算过程中，采用特定的定权方法，它可以固定所有不需要修正的常数元素，但协方差阵结构特殊不具普适性，在某些案例中不适用，而且不能解决重复元素的问题;Schaffrin［7，8］等提出多元整体最小二乘法，能避免系数矩阵出现重复元素，但是不能解决某些常数元素不需要修正的问题，且没有考虑观测值的权;Shen 等［9］提出用拉格朗日乘数法解算加权整体最小二乘，计算过程更简单，但他们仍然采用Schaffrin 等提出的定权方法，还是没有考虑系数矩阵元素的重复性;Mahboub［10］提出对系数矩阵的协方差阵加入特殊说明来自动获得精确解，但是该方法过程复杂，计算所需时间长，收敛速度较慢;袁庆等［6］将加权整体最小二乘法应用于三维基准转换中，达到改正系数矩阵所有数据元素和固定所有常数元素的作用，但是该方法没有考虑系数矩阵中的重复元素，相等元素的改正数不等。本文引用文献［9］的迭代方法且结合Mahboub 提出的定权理论，推导出更具普适性、收敛速度更快、自动解决系数矩阵元素的重复性的迭代算法，并将此方法应用于直线拟合和三维小角度基准转换验证该算法的效果。

2 改进的加权整体最小二乘算法

关于系数矩阵定权比较成熟的方法当属QA构造法，将QA分解为Q0∈Rm×m、Qx∈Rn×n［8］:

式中，Qy∈Rn×n、QA∈Rmn×mn分别为观测向量的协方差矩阵和系数矩阵列向量化后协方差矩阵，Py∈Rn×n、PA∈Rmn×mn分别为观测向量是权阵和系数矩阵列向量化后的权阵。它将PA看成是系数矩阵A列向量化后的向量权阵，使得PA可以对A 中的每一个元素定权，但是它却没有考虑元素之间的相关性，认为各个元素间是相互独立的。

实际应用发现，很多情况下系数矩阵中个元素间并不是完全相互独立的，某些元素重复出现，这些重复元素之间应该存在相关性，从理论的严密性来说文献［11］的方法不够严密，得出的解不是最优解。2012年Mahboub［10］提出利用一定原则构造考虑元素间相关性的协方差阵，该方法利用以下五条原则构造系数矩阵的协方差矩阵:

1)如果系数阵的某个元素重复出现，认为这两个元素100%相关，因此这两个元素之间的协方差等于重复元素的自方差;

2)假如系数阵的某个元素以其相反数的形式重复，认定这两个元素100%负相关，因此这两个元素之间的协方差等于重复元素的自方差的相反数;

3)如果系数阵的某个元素是常数，认为其方差为零;

4)系数阵中两个不同元素，若两者明显相关，他们的协方差即为其实际值，否则为零;

5)上述规则在同方差情况中同样适用，若元素是随机数只需用数字1 作为其方差，若是常数其自方差为零。

本文采用上述五条规则，重新构造系数矩阵的协方差阵，然后结合文献［9］的迭代方法提出改进的加权整体最小二乘法。用改进的加权整体最小二乘法解算方程组

式中ey，y∈Rn×m，A，EA∈Rn×m，ξ，δξ∈Rm×1。具体的解算步骤如下:

1)根据上述五条规则生成QA;

2)设置初始值EA(0)=0，ξ(0)=(ATA)－1ATy;

3)

其中eA为数矩阵列向量化后向量的随机误差，“⊗”为kronecker 积算子，如:

5)计算观测向量的残差矩阵ey和单位权中误差σ0:

3 改进的加权整体最小二乘法在直线拟合和三维空间基准变换中的应用

3.1 直线拟合

为了检验改进加权整体最小二乘算法的效果，引用文献［12］的实验数据，坐标观测值(Xi，Yi)及其相应的权(WXi，WYi)列于表1。直线拟合的EIV模型为式(11)，由式(11)可知系数矩阵的第一列是常数不需要修正，按上述五条规则构造系数矩阵的协方差矩阵QA为式(12)。

将表1 中的观测值分别用加权整体最小二乘法(WTLS)和本文方法(IWTLS)进行求解，参数估值列于表2。由表2 可以看出，为求得拟合参数的精确值，取阀值ε0=10－10，WTLS 法需要迭代14 次，而IWTLS 法只需7 次，说明本文方法可行且在本实验中收敛速度更快。

表1 坐标观测值及其相应权［12］Tab.1 Coordinate observations and their corresponding weights［12］

3.2 空间三维布尔沙转换模型

当旋转角是小角度或初值已知，且控制点数不小于3 个时，布尔沙转换模型可写成:

其中，(XS，YS，ZS)、(XT，YT，ZT)分别为原始坐标系和目标坐标系;(X0，Y0，Z0)为平移参数，εX、εY、εZ为旋转参数，μ 为尺度参数。由式(13)得出，系数矩阵中不仅有不需要改正的常数，还含有重复的随机元素。为了验证改进的加权整体最小二乘法能控制系数矩阵重复元素的改正数，我们将该方法应用到该模型，实验数据如表3。

分别用LS、LS-TLS、WTLS 以及IWTLS 法求出转换参数;分别用7 个公共点和全部公共点进行坐标转换求解参数及精度;最后将各方法求解的参数估值及精度列于表4 和表5;用7 个公共点求解坐标转换参数时系数矩阵的误差矩阵如表6。

表2 WTLS 法和IWTLS 法解算的参数估值Tab.2 Estimated parameters with the methods:WTLS and IWTLS

表3 制点坐标观测值及相应精度［11］(单位:m)Tab.3 Coordinate observations of control points and their corresponding accuracy［11］(unit:m)

表4 7 个点求解参数及精度比较(1、3、6、7、9、10、11 号点)Tab.4 Comparison among the calculated parameters and accuracies with piont 1，3，6，7，9，10，11

表5 全部公共点求解参数及精度比较Tab.5 Comparison among the calculated parameters and accuracies with all pionts

表6 IWTLS 法系数矩阵的残差阵EA(单位:m)Tab.6 Residuals matrix of coefficient matrix with IWTLS(unit:m)

当迭代阀值ε0取10－10时，文献［10］的方法需要迭代117 次，本文方法只需58 次，再次证明本方法收敛速度比文献［10］的方法快。由表4、5 可以得出，IWTLS 方法较LS 方法精度提高30%以上，比混合最小二乘法提高了20%以上，比原有的WTLS的精度提高了5%以上。由表6 可以看出在IWTLS中，系数矩阵的常数元素被固定，只修改了随机元素，系数矩阵残差阵中重复元素的改正数相等，与文献［11］中系数阵的残差阵相比更为严密，使得整体最小二乘法在解算EIV 模型时更为合理。

4 结束语

1)使用加权整体最小二乘法，引入系数矩阵及观测向量的权阵，一是可以顾及观测值精度不等的情况;二是可以对系数矩阵A 中的部分元素加以改正，使EIV 模型更加合理，使整体最小二乘法应用更为广泛。

2)引用最新的定权理论，对加权整体最小二乘法的系数阵的权阵加以改进，提出了改进的加权整体最小二乘法，并将该方法应用到直线拟合、三维小角度坐标转换中。计算精度较最小二乘法提高30%以上，较混合整体最小二乘法提高20%以上。IWTLS 方法考虑了系数矩阵中元素间的相关性，对系数矩阵中重复元素的改正数加以控制，使相同随机元素的改正数相等，在理论上更严密。至于系数矩阵中元素间的相关性是否可以忽略，在什么条件下忽略，不可忽略的话对结果的影响有多大还需要进一步研究。

3)参考丁克良等［13］提出的整体最小二乘应用条件:从参数估计的角度讲，只有在系数矩阵误差对最小奇异值的大小影响显著时，整体最小二乘法效果才比较明显;当系数矩阵的误差对最小奇异值的影响较小时，即系数矩阵比较精确时，就不适于采用整体最小二乘法。所以在实际应用中，我们要根据实际情况判断是否需要用整体最小二乘法。

1 Golub G H and van Loan.An analysis of the total leastsquares problem［J］.SIAM J Numer Anal.，1980，17:883－893.

2 魏木生.广义最小二乘问题的理论和计算［M］.北京:科学出版社，2006.(Wei Musheng.Generalized least squares theory and calculation［M］.Beijing:Science Press，2006)

3 孔建，姚宜斌，吴寒.整体最小二乘的迭代解法［J］.武汉大学学报(信息科学版)，2010，35(6):711－714.(Kong Jian，Yao Yibin and Wu Han.Iterative method for total least－squares［J］.Geomatics and Information Science of Wuhan University，2010，35(6):711－714)

4 鲁铁定，周世健.总体最小二乘的迭代解法［J］.武汉大学学报(信息科学版)，2010，35(11):1 351－1 354.(Lu Tieding and Zhou Shijian.An itreative algorithm for total least squares estimation［J］.Geomatics and Information Science of Wuhan University，2010，35(11):1 351－1 354)

5 陆珏，陈义，郑波.总体最小二乘方法在坐标转换中的应用［J］.大地测量与地球动力学，2008，(5):77－81.(Lu Jue，Chen Yi and Zheng Bo.Applying total least squares to three-dimensional datum transformation［J］.Journal of Geodesy and Geodynamics，2008，(5):77－81)

6 Schaffrin B and Wieser A.On weighted total least－squares adjustment for linear regression［J］.J Geodes.，2008，82(7):415 –421.

7 Schaffrin B and Felus Y A.On the multivariate total least squares approach to empirical coordinate transformation［J］.J Geodes.，2008，82(6):373－383.

8 Schaffrin B and Felus Y A.A multivariate total least-squares adjustment for empirical affine transformations［J］.International Association of Geodesy Symposia，2008，132(3):238－242.

9 Shen Y Z，Li B F and Chen Y.An iterative solution of weighted total least-squares adjustment［J］.Journal of Geodesy，2010，85(4):229－238.

10 Vahid Mahboub.On weighted total least－squares for geodetic transfor mations［J］.J Geod.，2012，86(5):359－367.

11 袁庆，楼立志，陈玮娴.加权总体最小二乘在三维基准中的应用［J］.测绘学报，2011，40:115－119.(Yuan Qing，Lou Lizhi and Chen Weixian.The application of the weighted total least-squares to three dimensional－datum transformation［J］.Acta Geodaetica et Cartographica Sinica，2011，40:115－119)

12 Neri F，Saitta G and Chiofalo S.Accurate and straightforward approach to line regression analysis of error-affected experimental data［J］.Phys E.，1989，22(4):215 –217.

13 丁克良，欧吉坤，陈义.整体最小二乘法及其在测量数据处理中的应用［A］.中国测绘学会第九次全国会员代表大会暨学会成立50 周年纪念大会论文集［C］.2009，399－405.(Ding Keliang，Ou Jikun and Chen Yi.Total leastsquares and its application in the measurement data processing［A］.Proceedings of the ninth national congress of the members of China surveying and mapping society and association established 50 anniversary conference［C］.2009，399－405)