楊仕平 范東明 龍玉春

1)西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院測(cè)繪工程系，成都 610031

2)重慶市合川龍市中學(xué)，合川401564

1 引言

從Golub 和van Loan［1］提出利用基于奇異值分解法(SVD，Singular Value Decom-position)的整體最小二乘法求解EIV 模型至今，整體最小二乘法一直在不斷地發(fā)展。鑒于SVD 方法計(jì)算復(fù)雜不利于編程，研究人員陸續(xù)提出了整體最小二乘的其他解法，如迭代算法、完全正交分解法等改進(jìn)方法［2－4］。

在大地測(cè)量和工程測(cè)量數(shù)據(jù)處理中，絕大部分都可以歸結(jié)為參數(shù)估計(jì)問題，這些參數(shù)估計(jì)問題大致可以分為普通最小二乘估計(jì)問題和整體最小二乘估計(jì)問題，而整體最小二乘估計(jì)問題可進(jìn)一步分為系數(shù)陣的元素全部含有隨機(jī)誤差及只有部分含有隨機(jī)誤差兩種情況。

針對(duì)測(cè)繪數(shù)據(jù)處理中系數(shù)矩陣有特殊結(jié)構(gòu)的EIV 模型，陸鈺等［5］提出用LS-TLS 將系數(shù)矩陣中不需要修正的列固定，但它不能固定所有的不需要修正的元素，不能考慮重復(fù)元素;Schaffrin 等［6］提出在加權(quán)整體最小二乘法解算過程中，采用特定的定權(quán)方法，它可以固定所有不需要修正的常數(shù)元素，但協(xié)方差陣結(jié)構(gòu)特殊不具普適性，在某些案例中不適用，而且不能解決重復(fù)元素的問題;Schaffrin［7，8］等提出多元整體最小二乘法，能避免系數(shù)矩陣出現(xiàn)重復(fù)元素，但是不能解決某些常數(shù)元素不需要修正的問題，且沒有考慮觀測(cè)值的權(quán);Shen 等［9］提出用拉格朗日乘數(shù)法解算加權(quán)整體最小二乘，計(jì)算過程更簡(jiǎn)單，但他們?nèi)匀徊捎肧chaffrin 等提出的定權(quán)方法，還是沒有考慮系數(shù)矩陣元素的重復(fù)性;Mahboub［10］提出對(duì)系數(shù)矩陣的協(xié)方差陣加入特殊說明來自動(dòng)獲得精確解，但是該方法過程復(fù)雜，計(jì)算所需時(shí)間長，收斂速度較慢;袁慶等［6］將加權(quán)整體最小二乘法應(yīng)用于三維基準(zhǔn)轉(zhuǎn)換中，達(dá)到改正系數(shù)矩陣所有數(shù)據(jù)元素和固定所有常數(shù)元素的作用，但是該方法沒有考慮系數(shù)矩陣中的重復(fù)元素，相等元素的改正數(shù)不等。本文引用文獻(xiàn)［9］的迭代方法且結(jié)合Mahboub 提出的定權(quán)理論，推導(dǎo)出更具普適性、收斂速度更快、自動(dòng)解決系數(shù)矩陣元素的重復(fù)性的迭代算法，并將此方法應(yīng)用于直線擬合和三維小角度基準(zhǔn)轉(zhuǎn)換驗(yàn)證該算法的效果。

2 改進(jìn)的加權(quán)整體最小二乘算法

關(guān)于系數(shù)矩陣定權(quán)比較成熟的方法當(dāng)屬Q(mào)A構(gòu)造法，將QA分解為Q0∈Rm×m、Qx∈Rn×n［8］:

式中，Qy∈Rn×n、QA∈Rmn×mn分別為觀測(cè)向量的協(xié)方差矩陣和系數(shù)矩陣列向量化后協(xié)方差矩陣，Py∈Rn×n、PA∈Rmn×mn分別為觀測(cè)向量是權(quán)陣和系數(shù)矩陣列向量化后的權(quán)陣。它將PA看成是系數(shù)矩陣A列向量化后的向量權(quán)陣，使得PA可以對(duì)A 中的每一個(gè)元素定權(quán)，但是它卻沒有考慮元素之間的相關(guān)性，認(rèn)為各個(gè)元素間是相互獨(dú)立的。

實(shí)際應(yīng)用發(fā)現(xiàn)，很多情況下系數(shù)矩陣中個(gè)元素間并不是完全相互獨(dú)立的，某些元素重復(fù)出現(xiàn)，這些重復(fù)元素之間應(yīng)該存在相關(guān)性，從理論的嚴(yán)密性來說文獻(xiàn)［11］的方法不夠嚴(yán)密，得出的解不是最優(yōu)解。2012年Mahboub［10］提出利用一定原則構(gòu)造考慮元素間相關(guān)性的協(xié)方差陣，該方法利用以下五條原則構(gòu)造系數(shù)矩陣的協(xié)方差矩陣:

1)如果系數(shù)陣的某個(gè)元素重復(fù)出現(xiàn)，認(rèn)為這兩個(gè)元素100%相關(guān)，因此這兩個(gè)元素之間的協(xié)方差等于重復(fù)元素的自方差;

2)假如系數(shù)陣的某個(gè)元素以其相反數(shù)的形式重復(fù)，認(rèn)定這兩個(gè)元素100%負(fù)相關(guān)，因此這兩個(gè)元素之間的協(xié)方差等于重復(fù)元素的自方差的相反數(shù);

3)如果系數(shù)陣的某個(gè)元素是常數(shù)，認(rèn)為其方差為零;

4)系數(shù)陣中兩個(gè)不同元素，若兩者明顯相關(guān)，他們的協(xié)方差即為其實(shí)際值，否則為零;

5)上述規(guī)則在同方差情況中同樣適用，若元素是隨機(jī)數(shù)只需用數(shù)字1 作為其方差，若是常數(shù)其自方差為零。

本文采用上述五條規(guī)則，重新構(gòu)造系數(shù)矩陣的協(xié)方差陣，然后結(jié)合文獻(xiàn)［9］的迭代方法提出改進(jìn)的加權(quán)整體最小二乘法。用改進(jìn)的加權(quán)整體最小二乘法解算方程組

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

1)根據(jù)上述五條規(guī)則生成QA;

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

3)

其中eA為數(shù)矩陣列向量化后向量的隨機(jī)誤差，“?”為kronecker 積算子，如:

5)計(jì)算觀測(cè)向量的殘差矩陣ey和單位權(quán)中誤差σ0:

3 改進(jìn)的加權(quán)整體最小二乘法在直線擬合和三維空間基準(zhǔn)變換中的應(yīng)用

3.1 直線擬合

為了檢驗(yàn)改進(jìn)加權(quán)整體最小二乘算法的效果，引用文獻(xiàn)［12］的實(shí)驗(yàn)數(shù)據(jù)，坐標(biāo)觀測(cè)值(Xi，Yi)及其相應(yīng)的權(quán)(WXi，WYi)列于表1。直線擬合的EIV模型為式(11)，由式(11)可知系數(shù)矩陣的第一列是常數(shù)不需要修正，按上述五條規(guī)則構(gòu)造系數(shù)矩陣的協(xié)方差矩陣QA為式(12)。

將表1 中的觀測(cè)值分別用加權(quán)整體最小二乘法(WTLS)和本文方法(IWTLS)進(jìn)行求解，參數(shù)估值列于表2。由表2 可以看出，為求得擬合參數(shù)的精確值，取閥值ε0=10－10，WTLS 法需要迭代14 次，而IWTLS 法只需7 次，說明本文方法可行且在本實(shí)驗(yàn)中收斂速度更快。

表1 坐標(biāo)觀測(cè)值及其相應(yīng)權(quán)［12］Tab.1 Coordinate observations and their corresponding weights［12］

3.2 空間三維布爾沙轉(zhuǎn)換模型

當(dāng)旋轉(zhuǎn)角是小角度或初值已知，且控制點(diǎn)數(shù)不小于3 個(gè)時(shí)，布爾沙轉(zhuǎn)換模型可寫成:

其中，(XS，YS，ZS)、(XT，YT，ZT)分別為原始坐標(biāo)系和目標(biāo)坐標(biāo)系;(X0，Y0，Z0)為平移參數(shù)，εX、εY、εZ為旋轉(zhuǎn)參數(shù)，μ 為尺度參數(shù)。由式(13)得出，系數(shù)矩陣中不僅有不需要改正的常數(shù)，還含有重復(fù)的隨機(jī)元素。為了驗(yàn)證改進(jìn)的加權(quán)整體最小二乘法能控制系數(shù)矩陣重復(fù)元素的改正數(shù)，我們將該方法應(yīng)用到該模型，實(shí)驗(yàn)數(shù)據(jù)如表3。

分別用LS、LS-TLS、WTLS 以及IWTLS 法求出轉(zhuǎn)換參數(shù);分別用7 個(gè)公共點(diǎn)和全部公共點(diǎn)進(jìn)行坐標(biāo)轉(zhuǎn)換求解參數(shù)及精度;最后將各方法求解的參數(shù)估值及精度列于表4 和表5;用7 個(gè)公共點(diǎn)求解坐標(biāo)轉(zhuǎn)換參數(shù)時(shí)系數(shù)矩陣的誤差矩陣如表6。

表2 WTLS 法和IWTLS 法解算的參數(shù)估值Tab.2 Estimated parameters with the methods:WTLS and IWTLS

表3 制點(diǎn)坐標(biāo)觀測(cè)值及相應(yīng)精度［11］(單位:m)Tab.3 Coordinate observations of control points and their corresponding accuracy［11］(unit:m)

表4 7 個(gè)點(diǎn)求解參數(shù)及精度比較(1、3、6、7、9、10、11 號(hào)點(diǎn))Tab.4 Comparison among the calculated parameters and accuracies with piont 1，3，6，7，9，10，11

表5 全部公共點(diǎn)求解參數(shù)及精度比較Tab.5 Comparison among the calculated parameters and accuracies with all pionts

表6 IWTLS 法系數(shù)矩陣的殘差陣EA(單位:m)Tab.6 Residuals matrix of coefficient matrix with IWTLS(unit:m)

當(dāng)?shù)y值ε0取10－10時(shí)，文獻(xiàn)［10］的方法需要迭代117 次，本文方法只需58 次，再次證明本方法收斂速度比文獻(xiàn)［10］的方法快。由表4、5 可以得出，IWTLS 方法較LS 方法精度提高30%以上，比混合最小二乘法提高了20%以上，比原有的WTLS的精度提高了5%以上。由表6 可以看出在IWTLS中，系數(shù)矩陣的常數(shù)元素被固定，只修改了隨機(jī)元素，系數(shù)矩陣殘差陣中重復(fù)元素的改正數(shù)相等，與文獻(xiàn)［11］中系數(shù)陣的殘差陣相比更為嚴(yán)密，使得整體最小二乘法在解算EIV 模型時(shí)更為合理。

4 結(jié)束語

1)使用加權(quán)整體最小二乘法，引入系數(shù)矩陣及觀測(cè)向量的權(quán)陣，一是可以顧及觀測(cè)值精度不等的情況;二是可以對(duì)系數(shù)矩陣A 中的部分元素加以改正，使EIV 模型更加合理，使整體最小二乘法應(yīng)用更為廣泛。

2)引用最新的定權(quán)理論，對(duì)加權(quán)整體最小二乘法的系數(shù)陣的權(quán)陣加以改進(jìn)，提出了改進(jìn)的加權(quán)整體最小二乘法，并將該方法應(yīng)用到直線擬合、三維小角度坐標(biāo)轉(zhuǎn)換中。計(jì)算精度較最小二乘法提高30%以上，較混合整體最小二乘法提高20%以上。IWTLS 方法考慮了系數(shù)矩陣中元素間的相關(guān)性，對(duì)系數(shù)矩陣中重復(fù)元素的改正數(shù)加以控制，使相同隨機(jī)元素的改正數(shù)相等，在理論上更嚴(yán)密。至于系數(shù)矩陣中元素間的相關(guān)性是否可以忽略，在什么條件下忽略，不可忽略的話對(duì)結(jié)果的影響有多大還需要進(jìn)一步研究。

3)參考丁克良等［13］提出的整體最小二乘應(yīng)用條件:從參數(shù)估計(jì)的角度講，只有在系數(shù)矩陣誤差對(duì)最小奇異值的大小影響顯著時(shí)，整體最小二乘法效果才比較明顯;當(dāng)系數(shù)矩陣的誤差對(duì)最小奇異值的影響較小時(shí)，即系數(shù)矩陣比較精確時(shí)，就不適于采用整體最小二乘法。所以在實(shí)際應(yīng)用中，我們要根據(jù)實(shí)際情況判斷是否需要用整體最小二乘法。

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

2 魏木生.廣義最小二乘問題的理論和計(jì)算［M］.北京:科學(xué)出版社，2006.(Wei Musheng.Generalized least squares theory and calculation［M］.Beijing:Science Press，2006)

3 孔建，姚宜斌，吳寒.整體最小二乘的迭代解法［J］.武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，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］.武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，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 陸玨，陳義，鄭波.總體最小二乘方法在坐標(biāo)轉(zhuǎn)換中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，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 袁慶，樓立志，陳瑋嫻.加權(quán)總體最小二乘在三維基準(zhǔn)中的應(yīng)用［J］.測(cè)繪學(xué)報(bào)，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 丁克良，歐吉坤，陳義.整體最小二乘法及其在測(cè)量數(shù)據(jù)處理中的應(yīng)用［A］.中國測(cè)繪學(xué)會(huì)第九次全國會(huì)員代表大會(huì)暨學(xué)會(huì)成立50 周年紀(jì)念大會(huì)論文集［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)