加權整體最小二乘算法的改進*

楊仕平 范東明 龍玉春

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)