楊勇喜,賈東振,何秀鳳

(河海大學(xué) 衛(wèi)星及空間信息應(yīng)用研究所，江蘇 南京 210098)

基于選權(quán)迭代法的抗差整體最小二乘及其應(yīng)用

楊勇喜,賈東振,何秀鳳

(河海大學(xué) 衛(wèi)星及空間信息應(yīng)用研究所，江蘇 南京 210098)

在測量數(shù)據(jù)處理中，觀測向量與系數(shù)矩陣同時(shí)存在偶然誤差時(shí)，整體最小二乘法能夠得到更高精度的參數(shù)解,但整體最小二乘法無抗差能力，觀測向量中的粗差將對參數(shù)求解產(chǎn)生較大影響。為解決上述問題，采用拉格朗日極值法推導(dǎo)了基于選權(quán)迭代法的抗差整體最小二乘計(jì)算公式,通過三維坐標(biāo)轉(zhuǎn)換參數(shù)求解實(shí)例對3種選權(quán)迭代法進(jìn)行分析。結(jié)果表明,IGG法在抗差整體最小二乘解法中抗差效果最好。

粗差；抗差整體最小二乘；選權(quán)迭代；三維坐標(biāo)轉(zhuǎn)換；IGG法

測量平差參數(shù)求解中，觀測向量僅存在偶然誤差時(shí)，可采用最小二乘法(LS)求得參數(shù)的最優(yōu)解。若觀測向量存在粗差，則采用抗差最小二乘法(抗差LS)剔除或減弱粗差的影響[1]。對于系數(shù)矩陣和觀測向量同時(shí)存在偶然誤差的情況，需要采用整體最小二乘方法。整體最小二乘問題首先由Golub和Van Loan[2]提出，并給出了奇異值分解的解法。后來Schaffrin提出了基于拉格朗日求極值的迭代法[3]。魏木生[4]、俞錦成[5]等對整體最小二乘理論進(jìn)行了深入細(xì)致的研究。實(shí)際應(yīng)用領(lǐng)域，如線性回歸[6]、空間后方交會(huì)[7]、建筑物沉降預(yù)測[8]等，整體最小二乘法均能夠得到更合理、更高精度的解。但是整體最小二乘方法沒有抗差能力，如果觀測向量存在粗差，參數(shù)估計(jì)將會(huì)偏離實(shí)際[9-10]。將選權(quán)迭代法應(yīng)用于整體最小二乘法的抗差估計(jì)，可以有效地抵御粗差對參數(shù)估計(jì)的影響。陳瑋嫻[11]、陳義[12]等都曾將選權(quán)迭代法應(yīng)用于三維坐標(biāo)轉(zhuǎn)換的整體最小二乘估計(jì)，但是選權(quán)迭代法有多種形式，哪種形式最有效，現(xiàn)有文獻(xiàn)鮮有提及。本文推導(dǎo)了基于選權(quán)迭代法的抗差整體最小二乘計(jì)算公式，并通過三維坐標(biāo)轉(zhuǎn)換對不同的選權(quán)迭代法進(jìn)行了對比和分析。結(jié)果表明，IGG法抗差效果最好。

1 數(shù)學(xué)模型

1.1 抗差整體最小二乘法

對于觀測方程L=AX，當(dāng)觀測向量L與系數(shù)矩陣A獨(dú)立時(shí)，EIV(Errors-In-Variables)模型[13]為

L-eL=(A-EA)X.

(1)

式中：L∈Rm×1為觀測向量，eL∈Rm×1為觀測向量的偶然誤差，A∈Rm×n為系數(shù)矩陣(列滿秩)，EA∈Rm×n為系數(shù)矩陣的偶然誤差，X∈Rn×1為待求參數(shù)。對應(yīng)的隨機(jī)模型

(2)

L-AX0-Aδx-eL+BeA=0.

(3)

式中B=(X0)T?Im，Im∈Rm×m。

整體最小二乘的目標(biāo)函數(shù)為

(4)

(5)

式中，ρ(eL(i))為M估計(jì)中選取的函數(shù)。利用拉格朗日求極值的方法，建立拉格朗日函數(shù)[11]

2λT(L-AX0-Aδx-eL+((X0)T?Im)eA).

(6)

對上式求eL偏導(dǎo)數(shù)，并令其為0。

(7)

(8)

類似地，式(6)分別對eA,λ,δx求偏導(dǎo)數(shù)，令其為0，并求出極值點(diǎn)

(9)

L-AX0-Aδx-eL+BeA=0,

(10)

(11)

(12)

(13)

由以上推導(dǎo)，結(jié)合文獻(xiàn)[11]，得到如下的計(jì)算步驟(上標(biāo)i為第i步迭代)：

1)設(shè)置初值

2)從i=0開始，依次計(jì)算

?Im,

(14)

根據(jù)式(13)得

(15)

并計(jì)算

X(i+1)=X(i)+δx(i+1),

(16)

(17)

(18)

1.2 常用選權(quán)迭代法

選權(quán)迭代法在抗差估計(jì)中應(yīng)用廣泛，常用的有以下3種形式。

1.2.1 胡貝爾法(Huber)

胡貝爾法確定的權(quán)因子為

(19)

當(dāng)所有改正數(shù)均在-c和c之間時(shí)，胡貝爾估計(jì)就是經(jīng)典的最小二乘估計(jì)。而當(dāng)改正數(shù)大于c時(shí)，其w(v)與改正數(shù)成反比，v越大，對應(yīng)的w(v)越小，權(quán)也越小，與此相應(yīng)該觀測值對參數(shù)估計(jì)的影響也越小。

1.2.2 漢佩爾法(Hampel)

漢佩爾法確定的權(quán)因子為

(20)

1.2.3 IGG法

IGG法是周江文在1989年提出的一種抗差權(quán)函數(shù)構(gòu)造方法[14]，對應(yīng)的權(quán)因子為

(21)

2 算例及分析

采用文獻(xiàn)[15]中的七參數(shù)坐標(biāo)轉(zhuǎn)換模型，求解WGS-84坐標(biāo)系和北京54坐標(biāo)系下三維坐標(biāo)的轉(zhuǎn)換參數(shù)。根據(jù)前面的求解步驟，計(jì)算分析不同選權(quán)迭代形式的抗差整體最小二乘法的結(jié)果。

表1給出了5個(gè)已知點(diǎn)在WGS-84和北京54坐標(biāo)系下的坐標(biāo)值，根據(jù)布爾莎模型求解WGS-84和北京54坐標(biāo)系之間的轉(zhuǎn)換參數(shù)。

首先采用最小二乘方法和整體最小二乘方法對七參數(shù)求解，然后在#1點(diǎn)的X坐標(biāo)加入3 m的粗差，再采用兩種方法求解，結(jié)果見表2。

表1 已知點(diǎn)在WGS-84坐標(biāo)系和北京54坐標(biāo)系下的坐標(biāo)值 m

表2 LS和TLS參數(shù)計(jì)算結(jié)果

由表2可知，未加入粗差時(shí)整體最小二乘方法的估計(jì)結(jié)果比最小二乘方法的精度高。加入粗差后，最小二乘方法和整體最小二乘方法均嚴(yán)重偏離了結(jié)果。

采用抗差整體最小二乘方法對加入粗差后的坐標(biāo)進(jìn)行參數(shù)估計(jì)?？共罟烙?jì)中的選權(quán)迭代方法依次采用胡貝爾法、漢佩爾法和IGG法。對于胡貝爾法設(shè)計(jì)4個(gè)方案。方案1：c=σ；方案2：c=1.5σ；方案3：c=2σ；方案4：c=2.5σ。參數(shù)計(jì)算結(jié)果如表3所示。

表3 胡貝爾法各方案比較

由表3可以看出，方案3即c=2σ時(shí)抗差估計(jì)效果最佳，當(dāng)c=2.5σ時(shí)與含粗差的整體最小二乘方法計(jì)算結(jié)果一致。

對于漢佩爾方法也設(shè)計(jì)4個(gè)方案。方案1：a=σ,b=1.5σ,c=3σ；方案2：a=1.5σ,b=2σ,c=4σ；方案3：a=2σ,b=2.5σ,c=5σ；方案4：a=2.5σ,b=3σ,c=6σ。參數(shù)計(jì)算結(jié)果如表4所示，方案1的抗差估計(jì)效果最佳。對于IGG法，k依次取1,10-1,10-2,10-4。參數(shù)估計(jì)結(jié)果如表5所示。

表4 漢佩爾法各方案比較

表5 IGG法各方案比較

由表5中數(shù)據(jù)可知，k的取值對結(jié)果影響很小，這與文獻(xiàn)[1]中的論述是一致的。

由上面的計(jì)算可知，3種選權(quán)迭代方法都有抵御粗差的能力，胡貝爾法方案3與漢佩爾法方案1的結(jié)果接近，但精度稍差；IGG法精度最高。IGG法從測量誤差理論來看，如果數(shù)據(jù)中只含有偶然誤差，界限mσ中m取1.5(在±1.5σ之外的概率僅為0.13)，這個(gè)區(qū)間以外的觀測值既不能完全排除又要限制其有害的影響。當(dāng)超過2.5σ(概率為0.01)，認(rèn)為是粗差，給予淘汰。IGG方案充分考慮了測量數(shù)據(jù)的實(shí)際情況，是一種十分有效的抗差方案。

3 結(jié) 論

為解決觀測向量中粗差對整體最小二乘結(jié)果產(chǎn)生顯著影響的問題，本文首先采用拉格朗日極值法推導(dǎo)了基于選權(quán)迭代法的抗差整體最小二乘計(jì)算公式，然后通過三維坐標(biāo)轉(zhuǎn)換參數(shù)求解實(shí)例對3種選權(quán)迭代法進(jìn)行了分析，得出以下結(jié)論：

1)在不含粗差情況下，整體最小二乘法比最小二乘法有更高精度的參數(shù)解，但是整體最小二乘法沒有抵御粗差的能力，需要加入抗差算法。

2)選權(quán)迭代法應(yīng)用于整體最小二乘法，可以較好地抵御粗差的影響，但是抗差能力有差異。計(jì)算結(jié)果表明，胡貝爾法與漢佩爾法抗差能力相當(dāng)，IGG法抗差能力最好。

[1]周江文,黃幼才,楊元喜. 抗差最小二乘法[M]. 武漢：華中理工大學(xué)出版社,1997.

[2]GOLUB G H,VAN LOAN C F. An Analysis of the Total Least Squares Problem[J]. SIAM Journal on Numerical Analysis.1980,17(6):883-893.

[3]SCHAFFRIM B. A note on Constrained Total Least-Squares Estimation[J]. Linear Algebra and Its Applications,2006,417:245-258.

[4]魏木生. 廣義最小二乘問題的理論和計(jì)算[M]. 北京:科學(xué)出版社,2006.

[5]俞錦成. 關(guān)于整體最小二乘問題的可解性[J]. 南京師范大學(xué)學(xué)報(bào):自然科學(xué)版,1996(1):12-16 .

[6]魯鐵定,陶本藻,周世健. 基于整體最小二乘法的線性回歸建模和解法[J]. 武漢大學(xué)學(xué)報(bào):信息科學(xué)版, 2008,33(5):504-507.

[7]陳義,陸玨,鄭波. 總體最小二乘方法在空間后方交會(huì)中的應(yīng)用[J]. 武漢大學(xué)學(xué)報(bào):信息科學(xué)版,2008,33(12):1271-1274.

[8]袁豹,岳東杰, 李成仁. 基于總體最小二乘的改進(jìn) GM (1, 1) 模型及其在建筑物沉降預(yù)測中的應(yīng)用[J]. 測繪工程, 2013, 22(3): 52-55.

[9]龔循強(qiáng)，劉國祥，李志林，等.總體最小二乘擬合問題求解方法的比較研究[J].測繪科學(xué)，2014，39(9)：29-33.

[10]鄧興升，孫虹虹，湯仲安.高斯牛頓迭代法解算非線性Bursa-Wolf模型的精度分析[J].測繪科學(xué)，2014，39(5)：93-95.

[11]陳瑋嫻,袁慶. 抗差總體最小二乘方法[J]. 大地測量與地球動(dòng)力學(xué),2012,32(6):111-114.

[12]陳義,陸玨. 以三維坐標(biāo)轉(zhuǎn)換為例解算穩(wěn)健總體最小二乘方法[J].測繪學(xué)報(bào),2012,41(5):715-722.

[13]YUNZHONG SHEN, BOFENG LI, YI CHEN. An iterative solution of weighted total least-squares adjustment[J]. Journal of Geodesy. 2011,85(4): 229-238.

[14]張勤,張菊清,岳東杰. 近代測量數(shù)據(jù)處理與應(yīng)用[M] 北京:測繪出版社,2011.

[15]武漢大學(xué)測繪學(xué)院測量平差學(xué)科組. 誤差理論與測量平差基礎(chǔ)[M]. 武漢:武漢大學(xué)出版社,2003.

[責(zé)任編輯：劉文霞]

Robust total least-squares based on selecting weight iteration method and its application

YANG Yong-xi，JIA Dong-zhen，HE Xiu-feng

(Institute of Satellite Navigation and Spatial Information System, Hohai University,Nanjing 210098,China)

In the measurement data processing,total least squares method can get more accurate parameter solution when observation vector and the coefficient matrix all exist random errors. But total least squares method doesn’t have the ability to resist the gross errors,which will have great impact on parameters. To solve the above problem,it minimizs the target function of classical Lagrange approach and deduces the robust total least-squares estimation calculation formula based on selecting weight iteration.Then the three selecting weight iteration methods are applied to three-dimensional coordinate transformation. The results show that IGG method works best in the robust total least-squares method.

gross errors; robust total least-squares method; selecting weight iteration; 3D coordinate transformation; IGG method

2013-11-27;補(bǔ)充更新日期：2014-10-15

國家自然科學(xué)基金資助項(xiàng)目(41274017)；江蘇省科技支撐計(jì)劃資助項(xiàng)目(BE2010316)

楊勇喜(1989-)，男，碩士研究生.

P207

：A

：1006-7949(2014)12-0056-04