陶葉青，高井祥，姚一飛

1. 淮海工學(xué)院測繪工程學(xué)院， 江蘇 連云港 222005； 2. 中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇 徐州 221008

?

基于中位數(shù)法的抗差總體最小二乘估計(jì)

陶葉青1，高井祥2，姚一飛2

1. 淮海工學(xué)院測繪工程學(xué)院， 江蘇 連云港 222005； 2. 中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇 徐州 221008

Foundation support： The National Natural Science Foundation of China (No. 41074010)；Public Science and Technology Research Funds Projects of Ocean(No. 201105004);Science Foundation of Huaihai Institute of Technology (No. Z2015006)

摘要：針對現(xiàn)有總體最小二乘抗差算法存在的缺陷，應(yīng)用中位數(shù)法確定模型參數(shù)的初值，提出了對模型的觀測向量與系數(shù)矩陣中的觀測元素進(jìn)行分類定權(quán)的思想，避免了中誤差估計(jì)偏差與隨機(jī)模型誤差對等價權(quán)函數(shù)抗差性的影響。基于中位數(shù)法建立總體最小二乘抗差迭代算法，并結(jié)合算例對算法進(jìn)行驗(yàn)證。結(jié)果表明，在相同觀測樣本條件下，本文提出的算法擬合的精度高于傳統(tǒng)算法擬合的精度。

關(guān)鍵詞：總體最小二乘；抗差估計(jì)；中位數(shù)

應(yīng)用總體最小二乘準(zhǔn)則(total least squares, TLS)能夠解決變量中含有誤差(error in variables, EIV)的模型參數(shù)估計(jì)問題，對其算法的研究在數(shù)據(jù)處理領(lǐng)域倍受關(guān)注[1-4]。近年來，趨于總體最小二乘算法的成熟，總體最小二乘理論在測繪工作中得到了廣泛應(yīng)用[5-8]。特別是能夠顧及模型隨機(jī)性質(zhì)算法的提出[9]，使得總體最小二乘模型在測繪數(shù)據(jù)處理中的應(yīng)用更加廣泛。同時，測繪工作者將經(jīng)典測量平差模型中廣泛存在的參數(shù)估計(jì)、精度評定、不適定模型的正則化、抗差估計(jì)、附有約束條件的模型參數(shù)估計(jì)等問題系統(tǒng)性地引入總體最小二乘準(zhǔn)則下進(jìn)行討論[10-17]：文獻(xiàn)[11]對總體最小二乘模型單位權(quán)中誤差的無偏估計(jì)進(jìn)行研究，并建立單位權(quán)中誤差無偏估計(jì)的迭代算法；文獻(xiàn)[13]對病態(tài)總體最小二模型的正則化算法進(jìn)行研究，建立了病態(tài)模型的廣義正則化算法；文獻(xiàn)[15]應(yīng)用高斯-赫爾墨特模型建立總體最小二乘的抗差估計(jì)算法；文獻(xiàn)[16]對附有等式約束的總體最小二乘算法進(jìn)行研究，并應(yīng)用極值函數(shù)解決參數(shù)估計(jì)問題，文獻(xiàn)[17]在此基礎(chǔ)上，應(yīng)用窮舉法建立附有不等式約束的總體最小二乘模型迭代算法。

現(xiàn)階段，總體最小二乘理論在測繪專業(yè)的應(yīng)用中取得豐富研究成果的同時，在抗差總體最小二乘理論方面缺乏適用性的研究成果，具有抗差性的總體最小二乘普遍算法是將經(jīng)典的等價權(quán)函數(shù)思想引入總體最小二乘理論[18]。直接將經(jīng)典的抗差理論應(yīng)用至總體最小二乘模型，筆者認(rèn)為存在兩個方面的缺陷：①EIV模型的系數(shù)矩陣中的元素與觀測向量中的元素精度不等，統(tǒng)一根據(jù)單位權(quán)中誤差確定它們的等價權(quán)函數(shù)的域值并不合理；②應(yīng)用總體最小二乘理論進(jìn)行參數(shù)估計(jì)時，隨機(jī)模型存在誤差，通過調(diào)節(jié)觀測值的權(quán)值來實(shí)現(xiàn)抗差性，會放大隨機(jī)模型誤差對算法的影響，算法的有效性值得商榷。

針對現(xiàn)有算法存在的不足，提出應(yīng)用中位數(shù)法建立抗差總體最小二乘估計(jì)，根據(jù)EIV模型的系數(shù)矩陣與觀測向量的殘差分類確定等價權(quán)函數(shù)的域值，避免隨機(jī)模型誤差與觀測元素精度差異等因素對抗差估計(jì)的影響。

1基于等價權(quán)函數(shù)的總體最小二乘抗差估計(jì)

1.1總體最小二乘模型

變量中含有誤差的參數(shù)估計(jì)模型為

v=(B+EB)x-l

(1)

式中，l為n×1階觀測向量；v為觀測向量中含有的n×1階誤差向量；B為模型的n×m階系數(shù)矩陣；EB為系數(shù)矩陣中含有的n×m階誤差矩陣；x為模型m×1階參數(shù)向量。令b=vec(EB)，vec(·)為矩陣的拉直符號，表示將矩陣按列拉直所得到的列向量。模型中誤差向量的隨機(jī)模型為

(2)

v=Bx-l

(3)

為實(shí)現(xiàn)對模型的觀測向量與系數(shù)矩陣中含有的誤差同時進(jìn)行最小化的約束，建立總體最小二乘約束準(zhǔn)則

(4)

當(dāng)Q、QB為單位矩陣時，式(4)退化為等權(quán)估計(jì)的總體最小二乘準(zhǔn)則。令

QB=Q0?Qb

(5)

vTPv+bT(P0?Pb)b=min

(6)

令PB=P0?Pb，則式(6)表示為

vTPv+bTPBb=min

(7)

1.2基于等價權(quán)函數(shù)的抗差估計(jì)

(8)

(9)

(10)

式(8)—式(10)中，其他變量的含義與參數(shù)的具體求解方法參見文獻(xiàn)[9]。根據(jù)觀測向量與系數(shù)矩陣中改正數(shù)的估值，應(yīng)用權(quán)函數(shù)對其進(jìn)行重新定權(quán)，以IGG權(quán)函數(shù)為例[19]

(11)

上述算法是現(xiàn)有研究成果中，實(shí)現(xiàn)總體最小二乘抗差估計(jì)的普遍思路，其優(yōu)點(diǎn)是理論基礎(chǔ)成熟、算法的抗差性容易實(shí)現(xiàn)。但是結(jié)合總體最小二乘模型的特點(diǎn)，根據(jù)上文的分析，直接將經(jīng)典測量平差抗差理論引用至總體最小二乘平差，其缺點(diǎn)也是顯著的。為克服觀測向量與系數(shù)矩陣中的元素精度不等，以及隨機(jī)模型誤差對總體最小二乘抗差性的影響，基于中位數(shù)建立總體最小二乘抗差估計(jì)算法。

2基于中位數(shù)法的抗差總體最小二乘估計(jì)

(12)

(13)

(14)

(15)

(16)

應(yīng)用中位數(shù)法確定參數(shù)的估值，可以有效降低含有粗差的觀測值對參數(shù)估值的影響，相關(guān)的研究成果表明，其崩潰污染率可以達(dá)到50%[21]。根據(jù)模型的觀測向量與系數(shù)矩陣中觀測元素的改正數(shù)，應(yīng)用中位數(shù)法分類確定它們的中誤差，并進(jìn)行分類定權(quán)的總體最小二乘抗差算法，可以有效避免由于觀測向量與系數(shù)矩陣中的觀測元素實(shí)際精度不等，以及隨機(jī)模型誤差對權(quán)函數(shù)域值確定的影響。

3算例與分析

應(yīng)用直線方程擬合對論文提出的算法進(jìn)行驗(yàn)證，并與現(xiàn)有算法進(jìn)行比較。設(shè)直線方程為y=a0x+a1，方程參數(shù)的模擬值a0=0.5、a1=10。直線上的點(diǎn)共有10個，其x軸的坐標(biāo)值列于表1。根據(jù)點(diǎn)的x軸坐標(biāo)與方程參數(shù)的模擬值計(jì)算點(diǎn)的y軸坐標(biāo)，列于表1。點(diǎn)坐標(biāo)嚴(yán)格滿足模擬參數(shù)所確定的直線方程，可以將其作為觀測元素的真值。

表1　點(diǎn)坐標(biāo)的模擬值

假設(shè)有n個觀測值時，直線的觀測方程表示為

(17)

將1號點(diǎn)的x軸坐標(biāo)、2號點(diǎn)的y軸坐標(biāo)混入4～5 dm的粗差，點(diǎn)1、2、3、4的其余坐標(biāo)混入1～6 cm的隨機(jī)誤差，應(yīng)用混入誤差的點(diǎn)1、2、3、4作為觀測值求解模型參數(shù)，混入誤差的點(diǎn)坐標(biāo)列于表2。

表2　混入誤差的點(diǎn)坐標(biāo)

在混入的誤差中，x軸方向的誤差值高于y軸方向的誤差，以模擬總體最小二乘模型的觀測向量與系數(shù)矩陣中的觀測元素精度不等。其余點(diǎn)不混入誤差，作為檢核點(diǎn)對模型的擬合精度進(jìn)行檢核。

應(yīng)用提出的中位數(shù)法求解模型參數(shù)，將選取的觀測點(diǎn)1、2、3、4分為4組，即(1、2、3)、(1、2、4)、(1、3、4)、(2、3、4)，應(yīng)用拉格朗日迭代算法分別求解模型參數(shù)的初值。最終求解的基于中位數(shù)的估值列于表3。根據(jù)傳統(tǒng)的總體最小二乘抗差算法，應(yīng)用點(diǎn)1、2、3、4求解的模型參數(shù)列于表3。模型的觀測向量與系數(shù)矩陣的權(quán)矩陣的初始值均取為單位1，如以4個點(diǎn)同時求解為例，式(8)中協(xié)因數(shù)陣的取值分別為

根據(jù)表2中觀測點(diǎn)混入誤差的特征，此時定義的隨機(jī)模型存在誤差。分別根據(jù)拉格朗日算法與中位數(shù)法，經(jīng)過三次迭代計(jì)算獲得模型參數(shù)的估值、改正數(shù)的估值、與中誤差的估值列于表3。

表3　不同算法獲得的估值

注： 1表示拉格朗日算法，2表示中位數(shù)算法；系數(shù)矩陣中的常向量改正數(shù)為0，故只列出其觀測值組成的列向量對應(yīng)的改正數(shù)估值。中誤差0.146/0.074分別表示觀測向量與系數(shù)矩陣對應(yīng)的中誤差。

由于總體最小二乘準(zhǔn)則或者最小二乘準(zhǔn)則會對觀測值中的誤差進(jìn)行“平攤”，并且應(yīng)用拉格朗日算法得到的中誤差是有偏的[9]，表3的結(jié)果表明，如果根據(jù)傳統(tǒng)算法，應(yīng)用權(quán)函數(shù)并不能有效地根據(jù)含有粗差的觀測值對參數(shù)進(jìn)行抗差估計(jì)。應(yīng)用表3中由中位數(shù)法得到的各類估值，按照本文提出的方法，根據(jù)權(quán)函數(shù)對觀測向量與系數(shù)矩陣中的觀測元素進(jìn)行分類定權(quán)，繼續(xù)進(jìn)行迭代計(jì)算，并以最小范數(shù)解作為模型參數(shù)的最終估值(式(16))。

當(dāng)觀測樣本有限時，總體最小二乘是一種有偏估計(jì)[9]，因此，不以中誤差為精度的評定標(biāo)準(zhǔn)。根據(jù)總體最小二乘的幾何意義[22]，以檢核點(diǎn)(表1，點(diǎn)5、6、7、8、9、10)到擬合直線垂直距離的真值與擬合值之差ΔD，比較不同算法的擬合精度

(18)

式中，Δy為檢核點(diǎn)y坐標(biāo)的真值與擬合值的差值。不同算法擬合的精度列于表4。結(jié)果表明，論文提出的基于中位數(shù)的總體最小二乘抗差算法擬合精度高于傳統(tǒng)總體最小二乘抗差算法，提出的應(yīng)用中位數(shù)法對模型的觀測向量與系數(shù)矩陣中的觀測元素進(jìn)行分類定權(quán)的思想是可行的。在相同觀測樣本條件下，本文提出的抗差算法擬合的模型精度更高。

在上述算例中，為了使得應(yīng)用中位數(shù)法求解模型參數(shù)時，每組觀測值中都含有粗差，分別在不同的點(diǎn)中加入粗差。當(dāng)應(yīng)用中位數(shù)法求解模型參數(shù)時，每組觀測樣本數(shù)小于傳統(tǒng)算法，使得樣本的粗差污染率變大。如果只在上述觀測點(diǎn)中的一個觀測點(diǎn)的縱橫坐標(biāo)中加入粗差，論文提出的算法在擬合精度上更有優(yōu)勢。

表4　不同算法擬合的精度

注：1表示拉格朗日算法，2表示中位數(shù)算法。

4結(jié)論

現(xiàn)有總體最小二乘抗差算法無法顧及中誤差估計(jì)偏差與隨機(jī)模型誤差等因素對等價權(quán)函數(shù)的影響。為克服這一缺陷，筆者應(yīng)用中位數(shù)建立總體最小二乘的抗差迭代算法，提出對模型觀測向量與系數(shù)矩陣中的觀測值進(jìn)行分類定權(quán)的思想，可提高等價權(quán)函數(shù)的有效性。試驗(yàn)結(jié)果表明，在相同觀測樣本、粗差污染率較高的條件下，本文所提算法得到的擬合精度高于傳統(tǒng)算法。

參考文獻(xiàn)：

[1]PEARSON K. On Lines and Planes of Closest Fit to Systems of Points in Space[J]. Philosophical Magazine Series 6, 1901, 2(11): 559-572.

[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]VAN HUFFEL S, VANDEWALLE J. Analysis and Solution of the Nongeneric Total Least Squares Problem[J]. SIAM Journal on Matrix Analysis and Applications, 1988, 9(3): 360-372.

[4]SCHAFFRIN B, FELUS Y A. On Total Least-squares Adjustment with Constraints[J]. International Association of Geodesy Symposia, 2005, 128: 417-421.

[5]SCHAFFRIN B, LEE I, CHOI Y, et al. Total Least Squares (TLS) for Geodetic Straight-line and Plane Adjustment[J]. Bollettino di Geodesia e Scienze Affini, 2006, 65(3): 141-168.

[6]MAHBOUB V. On Weighted Total Least-squares for Geodetic Transformations[J]. Journal of Geodesy, 2012, 86(5): 359-367.

[7]SCHAFFRIN B, FELUS Y A. On the Multivariate Total Least-squares Approach to Empirical Coordinate Transformations: Three Algorithms[J]. Journal of Geodesy, 2008, 82(6): 373-383.

[8]KUPFERER S. Application of the Total Least-squares Technology with Geodetic Problem Definitions[D]. Karlsruhe: University of Karlsruhe, 2005.

[9]SCHAFFRIN B, WIESER A. On Weighted Total Least-squares Adjustment for Linear Regression[J]. Journal of Geodesy, 2008, 82(7): 415-421.

[10]FULLER W A. Properties of Some Estimators for the Errors-in-Variables Model[J]. The Annals of Statistics, 1980, 8(2): 407- 422.

[11]TONG Xiaohua, JIN Yanmin, LI Lingyun. An Improved Weighted Total Least Squares Method with Applications in Linear Fitting and Coordinate Transformation[J]. Journal of Surveying Engineering, 2011, 137(4): 120-128.

[12]SHEN Yunzhong, LI Bofeng, CHEN Yi. An Iterative Solution of Weighted Total Least-squares Adjustment[J]. Journal of Geodesy, 2011, 85(4): 229-238.

[13]葛旭明, 伍吉倉. 病態(tài)總體最小二乘問題的廣義正則化[J]. 測繪學(xué)報, 2012, 41(3): 372-377.

GE Xuming, WU Jicang. Generalized Regularization to III-posed Total Least Squares Problem[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(3): 372-377.

[14]SCHAFFRIN B, SNOW K. Total Least-squares Regularization of Tykhonov Type and an Ancient Racetrack in Corinth[J]. Linear Algebra and Its Applications, 2010, 432(8): 2061-2076.

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

CHEN Yi,LU Jue.Performing 3D Similarity Transformation by Robust Total Least Squares[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(5): 715-722.

[16]SCHAFFRIN B. A Note on Constrained Total Least Squares Estimation[J]. Linear Algebra and Its Applications, 2006, 417(1): 245-258.

[17]ZHANG Songlin, TONG Xiaohua, ZHANG Kun. A Solution to EIV Model with Inequality Constraints and Its Geodetic Applications[J]. Journal of Geodesy, 2013, 87(1): 23-28.

[18]TAO Yeqing, GAO Jinxiang, YAO Yifei. TLS Algorithm for GPS Height Fitting Based on Robust Estimation[J] Survey Review, 2014, 46(336): 184-188.

[19]周江文. 經(jīng)典誤差理論與抗差估計(jì)[J]. 測繪學(xué)報, 1989, 18(2): 115-120.

ZHOU Jiangwen. Classic Theory of Errors and Robust Estimation[J]. Acta Geodaetica et Cartographica Sinica, 1989, 18(2): 115-120.

[20]ROUSSEEUW P J, WAGNER J. Robust Regression with a Distributed Intercept Using Least Median of Squares[J]. Computational Statistics & Data Analysis, 1994, 17(1): 65-76.

[21]YANG Yuanxi. Robust Estimation of Geodetic Datum Transformation[J]. Journal of Geodesy, 1999, 73(5): 268-274.

[22]劉經(jīng)南, 曾文憲, 徐培亮. 整體最小二乘估計(jì)的研究進(jìn)展[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2013, 38(5): 505-512.

LIU Jingnan, ZENG Wenxian, XU Peiliang. Overview of Total Least Squares Methods[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5): 505-512.

(責(zé)任編輯：宋啟凡)

Solution for Robust Total Least Squares Estimation Based on Median Method

TAO Yeqing1，GAO Jingxiang2，YAO Yifei2

1. Department of Geomatic Engineering, Huaihai Institute of Technology, Lianyungang 222005, China; 2. School of Environment Science and Spatial Informatics, China University of Mining & Technology, Xuzhou 221008, China

Abstract：Because the present algorithms of total least squares for robust estimation have disadvantages, solution for computation primary model parameters based on median method is proposed. And to get over the influence that estimation error of random model and error of mean square has, computation weight matrix of observation vector and coefficient matrix separately are proposed. Iterative algorithm of robust total least squares is established based on median method, and to prove the proposed solution to be feasible, an instance is cited. The numerical results of the instance clearly demonstrate that the presented solution is more accurate than the traditional method for line fitting.

Key words：total least squares; robust estimation; median method

基金項(xiàng)目：國家自然科學(xué)基金(41074010)；國家海洋公益專項(xiàng)(201105004)；淮海工學(xué)院科研基金(Z2015006)

中圖分類號：P207

文獻(xiàn)標(biāo)識碼：A

文章編號：1001-1595(2016)03-0297-05

作者簡介：第一 陶葉青(1984—)，男，博士，講師，研究方向?yàn)镚NSS導(dǎo)航與測量數(shù)據(jù)處理。E-mail： yenneytao@163.com

收稿日期：2015-05-05

引文格式：陶葉青，高井祥，姚一飛.基于中位數(shù)法的抗差總體最小二乘估計(jì)[J].測繪學(xué)報，2016,45(3)：297-301. DOI:10.11947/j.AGCS.2016.20150234.

TAO Yeqing，GAO Jingxiang，YAO Yifei.Solution for Robust Total Least Squares Estimation Based on Median Method[J]. Acta Geodaetica et Cartographica Sinica,2016,45(3):297-301. DOI:10.11947/j.AGCS.2016.20150234.

修回日期： 2015-11-19

First author： TAO Yeqing(1984—), male, PhD, lecturer, majors in GNSS navigation and survey data processing.