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

陶葉青，高井祥，姚一飛

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

?

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

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

1. 淮海工學(xué)院測(cè)繪工程學(xué)院， 江蘇 連云港 222005； 2. 中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(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)

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

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

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

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

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

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

1.1總體最小二乘模型

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

v=(B+EB)x-l

(1)

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

(2)

v=Bx-l

(3)

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

(4)

當(dāng)Q、QB為單位矩陣時(shí)，式(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基于等價(jià)權(quán)函數(shù)的抗差估計(jì)

(8)

(9)

(10)

式(8)—式(10)中，其他變量的含義與參數(shù)的具體求解方法參見文獻(xiàn)[9]。根據(jù)觀測(cè)向量與系數(shù)矩陣中改正數(shù)的估值，應(yīng)用權(quán)函數(shù)對(duì)其進(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)典測(cè)量平差抗差理論引用至總體最小二乘平差，其缺點(diǎn)也是顯著的。為克服觀測(cè)向量與系數(shù)矩陣中的元素精度不等，以及隨機(jī)模型誤差對(duì)總體最小二乘抗差性的影響，基于中位數(shù)建立總體最小二乘抗差估計(jì)算法。

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

(12)

(13)

(14)

(15)

(16)

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

3算例與分析

應(yīng)用直線方程擬合對(duì)論文提出的算法進(jìn)行驗(yàn)證，并與現(xiàn)有算法進(jìn)行比較。設(shè)直線方程為y=a0x+a1，方程參數(shù)的模擬值a0=0.5、a1=10。直線上的點(diǎn)共有10個(gè)，其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ù)所確定的直線方程，可以將其作為觀測(cè)元素的真值。

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

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

(17)

將1號(hào)點(diǎn)的x軸坐標(biāo)、2號(hào)點(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作為觀測(cè)值求解模型參數(shù)，混入誤差的點(diǎn)坐標(biāo)列于表2。

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

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

應(yīng)用提出的中位數(shù)法求解模型參數(shù)，將選取的觀測(cè)點(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。模型的觀測(cè)向量與系數(shù)矩陣的權(quán)矩陣的初始值均取為單位1，如以4個(gè)點(diǎn)同時(shí)求解為例，式(8)中協(xié)因數(shù)陣的取值分別為

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

表3　不同算法獲得的估值

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

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

當(dāng)觀測(cè)樣本有限時(shí)，總體最小二乘是一種有偏估計(jì)[9]，因此，不以中誤差為精度的評(píng)定標(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ù)法對(duì)模型的觀測(cè)向量與系數(shù)矩陣中的觀測(cè)元素進(jìn)行分類定權(quán)的思想是可行的。在相同觀測(cè)樣本條件下，本文提出的抗差算法擬合的模型精度更高。

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

表4　不同算法擬合的精度

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

4結(jié)論

現(xiàn)有總體最小二乘抗差算法無法顧及中誤差估計(jì)偏差與隨機(jī)模型誤差等因素對(duì)等價(jià)權(quán)函數(shù)的影響。為克服這一缺陷，筆者應(yīng)用中位數(shù)建立總體最小二乘的抗差迭代算法，提出對(duì)模型觀測(cè)向量與系數(shù)矩陣中的觀測(cè)值進(jìn)行分類定權(quán)的思想，可提高等價(jià)權(quán)函數(shù)的有效性。試驗(yàn)結(jié)果表明，在相同觀測(cè)樣本、粗差污染率較高的條件下，本文所提算法得到的擬合精度高于傳統(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]葛旭明, 伍吉倉(cāng). 病態(tài)總體最小二乘問題的廣義正則化[J]. 測(cè)繪學(xué)報(bào), 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]. 測(cè)繪學(xué)報(bào), 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]. 測(cè)繪學(xué)報(bào), 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é)報(bào)(信息科學(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)目：國(guó)家自然科學(xué)基金(41074010)；國(guó)家海洋公益專項(xiàng)(201105004)；淮海工學(xué)院科研基金(Z2015006)

中圖分類號(hào)：P207

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

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

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

收稿日期：2015-05-05

引文格式：陶葉青，高井祥，姚一飛.基于中位數(shù)法的抗差總體最小二乘估計(jì)[J].測(cè)繪學(xué)報(bào)，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.