楚 彬，范東明，劉 波，秦 寧

(西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 610031)

基于EIV模型的穩(wěn)健估計(jì)

楚 彬，范東明，劉 波，秦 寧

(西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 610031)

EIV(error-in-variables)模型同時(shí)考慮觀測(cè)向量和系數(shù)矩陣的誤差，自提出以來便得到廣泛應(yīng)用。目前針對(duì)EIV模型的整體最小二乘解法(TLS)假設(shè)觀測(cè)值僅含有偶然誤差，當(dāng)觀測(cè)值存在粗差時(shí)其解并不是最優(yōu)的。文中通過選定合適的權(quán)函數(shù)，結(jié)合加權(quán)整體最小二乘迭代算法，導(dǎo)出基于EIV模型的穩(wěn)健整體最小二乘迭代解法(RTLS)。線性擬合實(shí)驗(yàn)表明，文中方法能對(duì)粗差進(jìn)行定位，且估計(jì)量受粗差影響較小，具有穩(wěn)健性。

EIV模型；整體最小二乘估計(jì)(TLS)；穩(wěn)健估計(jì)；穩(wěn)健整體最小二乘(RTLS)；線性擬合

經(jīng)典的Gauss-Markov模型假定函數(shù)模型已知、非隨機(jī)，并且認(rèn)為系數(shù)矩陣是可以精確求定的，僅假定觀測(cè)值向量包含隨機(jī)誤差[1]。在許多實(shí)際問題中如數(shù)字地面模型擬合、大地測(cè)量反演、GIS空間數(shù)據(jù)分析、滑坡監(jiān)測(cè)和坐標(biāo)變換等數(shù)學(xué)模型中，觀測(cè)向量和描述函數(shù)模型的系數(shù)矩陣均由觀測(cè)數(shù)據(jù)組成，兩者都包含隨機(jī)誤差[2]。這類平差模型稱為EIV(error-in-variables)模型。由于EIV模型同時(shí)考慮觀測(cè)向量和系數(shù)矩陣的誤差，因此采用經(jīng)典最小二乘估計(jì)方法(LS)對(duì)其進(jìn)行估計(jì)則是有偏的[3]。針對(duì)此模型，Golub等人于1980年提出了著名的奇異值分解算法(SVD)，并將其命名為整體最小二乘法(TLS)[4]，自此，整體最小二乘估計(jì)引起了各領(lǐng)域的廣泛關(guān)注。除奇異值分解算法外，還有經(jīng)典的拉格朗日方法[5]，也可作為非線性Gauss-Helmert模型處理[6]。以上的解法都是在觀測(cè)向量和系數(shù)矩陣僅存在隨機(jī)誤差的前提下進(jìn)行的，當(dāng)觀測(cè)向量和系數(shù)矩陣還存在粗差時(shí)，則模型歪曲，造成參數(shù)估計(jì)嚴(yán)重失實(shí)。

穩(wěn)健估計(jì)法能夠保證所估的參數(shù)不受或少受模型誤差(首先指的是粗差)的影響，主要用來發(fā)現(xiàn)粗差和對(duì)粗差進(jìn)行定位[7]。Choi等最早對(duì)EIV模型進(jìn)行穩(wěn)健整體最小二乘估計(jì)[8]，但是其方法基于傳統(tǒng)的假設(shè)檢驗(yàn)，僅適用于觀測(cè)向量或系數(shù)矩陣只含有一個(gè)粗差的情況。陳義等提出基于非線性Gauss-Helmert模型的選權(quán)迭代總體最小二乘解法[9]，雖能獲得穩(wěn)健的參數(shù)估計(jì)，但是其模型協(xié)方差陣結(jié)構(gòu)特殊，未考慮觀測(cè)量之間的相關(guān)性，在某些案例中不適用，而且不能解決重復(fù)元素的問題。本文通過選擇合適的權(quán)函數(shù)，結(jié)合Mahboub提出的改進(jìn)的加權(quán)整體最小二乘迭代解法(WTLS)[10]，導(dǎo)出穩(wěn)健整體最小二乘迭代解法(RTLS)。通過迭代，含粗差觀測(cè)值的權(quán)函數(shù)元素的值會(huì)逐步趨近于0，不含粗差觀測(cè)值的權(quán)函數(shù)元素變化不大，函數(shù)模型具有抗拒粗差的本領(lǐng)，待估參數(shù)受粗差的影響較小。

1 EIV模型和加權(quán)整體最小二乘解(WTLS)

1.1 EIV模型

觀測(cè)變量含有誤差的EIV模型的線性函數(shù)關(guān)系式為[11]

y-ey=(A-EA)x.

(1)

式中：y為含有隨機(jī)誤差ey的m維效應(yīng)量觀測(cè)向量，A為含有隨機(jī)誤差EA的m×n維自變量觀測(cè)值或其函數(shù)，x為n維待估參數(shù)向量。

隨機(jī)誤差具有如下統(tǒng)計(jì)性質(zhì)：

(2)

實(shí)際應(yīng)用發(fā)現(xiàn)，很多情況下系數(shù)矩陣中各元素間并不是完全相互獨(dú)立的，系數(shù)矩陣中的元素可能重復(fù)出現(xiàn)，或某個(gè)元素是以某些元素為自變量的函數(shù)。因此，本文采用Mahboub提出利用一定原則構(gòu)造系數(shù)矩陣的協(xié)方差矩陣[11]。

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)為其方差為0。

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

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

1.2 加權(quán)整體最小二乘解(WTLS)

WTLS估計(jì)準(zhǔn)則為

(3)

加權(quán)整體最小二乘實(shí)質(zhì)就是在極值條件式(3)和誤差方程式(1)的約束條件下求得X的最佳估值。按Lagrange乘數(shù)法求解，構(gòu)成目標(biāo)函數(shù)為[5]

2λ[y-Ax-ey+(xT?Im)eA].

(4)

式中：λ為n×1維拉格朗日乘數(shù)向量，Im為m×m維單位矩陣。

為求得Φ的極小值，將上式對(duì)ey，eA，λ和x求偏導(dǎo)，并令其為0。

(5)

(6)

(7)

(8)

由式(5)和式(6)可求得殘差預(yù)測(cè)值為

(9)

(10)

在導(dǎo)出殘差預(yù)測(cè)公式后，采用Mahboub提出的改進(jìn)的加權(quán)整體最小二乘解法求得待估參數(shù)x，其迭代過程如下：

1) 根據(jù)上文提到的五原則構(gòu)造系數(shù)矩陣的協(xié)方差矩陣QA。

2 穩(wěn)健整體最小二乘解(RTLS)

2.1 權(quán)函數(shù)的選取

1）實(shí)驗(yàn)的進(jìn)行不能顧此失彼，對(duì)于單樁豎向抗壓承載力實(shí)驗(yàn)來說，為了避免荷載過重造成儀器損壞，影響實(shí)驗(yàn)數(shù)據(jù)的準(zhǔn)確，必須進(jìn)行必要的操作過程設(shè)計(jì)規(guī)范和細(xì)節(jié)控制說明，比如：當(dāng)?shù)鼗車寥缆∑饡r(shí)，其觀測(cè)點(diǎn)的荷載-沉降曲線必然出現(xiàn)陡降，實(shí)驗(yàn)中要注意記錄陡降前后的荷載數(shù)值。此外，當(dāng)?shù)鼗A(chǔ)樁后一次的沉降幅度超過前一次沉降幅度的2倍時(shí)，表明沉降過程尚未穩(wěn)定，因此，要注意記錄前一級(jí)的荷載數(shù)值。最后，在細(xì)節(jié)把控上，荷載-沉降曲線一定要按照信號(hào)傳感器反饋的時(shí)間曲線進(jìn)行繪制，必要的情況下，增加輔助曲線，增加檢測(cè)結(jié)果的準(zhǔn)確性。

以上的WTLS解僅適用于觀測(cè)向量和系數(shù)矩陣僅含有偶然誤差的情況，當(dāng)其還含有粗差時(shí)則會(huì)引起模型歪曲，造成參數(shù)估計(jì)嚴(yán)重失實(shí)。為了達(dá)到估值穩(wěn)健的目的，在目標(biāo)函數(shù)式中應(yīng)加入合適的權(quán)函數(shù)，使它具有抗拒粗差的本領(lǐng)，而不受粗差的影響。權(quán)函數(shù)有多種，目前最為著名的有最小范數(shù)法、Huber權(quán)函數(shù)、Hampel權(quán)函數(shù)和Krarup(丹麥法)權(quán)函數(shù)等。本文以Huber權(quán)函數(shù)為例，并根據(jù)文獻(xiàn)[12]的建議對(duì)其稍作改進(jìn)，所得權(quán)函數(shù)如下：

(11)

2.2 穩(wěn)健整體最小二乘解(RTLS)

定理1：假設(shè)A為n×n維實(shí)對(duì)稱矩陣，那么則有n×n維正交矩陣S和對(duì)角矩陣U(U中對(duì)角線上的元素為A的特征值)，使得STAS=U,STS=In。

將QA進(jìn)行Schur分解可得QA=SUST。結(jié)合1.2中加權(quán)整體最小二乘解(WTLS)和文獻(xiàn)[13]的最小二乘穩(wěn)健估計(jì)方法，可得穩(wěn)健整體最小二乘解(RTLS)，迭代過程如下：

其中：Wy(i)和WA(i)由式(11)確定，(·)-0.5為矩陣·的算術(shù)平方根逆根。在迭代過程中，根據(jù)式(9)和式(10)可求得每次迭代的殘差預(yù)測(cè)值，然后將其帶入式(11)對(duì)權(quán)函數(shù)進(jìn)行更正。通過迭代，含粗差觀測(cè)值的權(quán)函數(shù)元素的值會(huì)逐步趨近于0，不含粗差觀測(cè)值的權(quán)函數(shù)元素變化不大。因此，此方法不僅可以對(duì)粗差進(jìn)行定位，而且所估參數(shù)受粗差影響較小，具有穩(wěn)健性。

3 線性擬合實(shí)驗(yàn)分析

3.1 觀測(cè)值僅含有偶然誤差

為了考察穩(wěn)健整體最小二乘方法(RTLS)的效果，通過編寫本文所討論方法的程序，對(duì)平面回歸算例進(jìn)行計(jì)算并與LS和WTLS方法進(jìn)行比較。算例如下：設(shè)有平面方程z=1.5+1.6x-1.2y，取滿足平面方程的8組數(shù)值(不含誤差)統(tǒng)計(jì)于表1，然后對(duì)x和y添加ε1∈[-0.3～0.3]的隨機(jī)誤差得到組成系數(shù)矩陣的觀測(cè)值，對(duì)z添加ε2∈[-0.2～0.2]的隨機(jī)誤差構(gòu)成與系數(shù)矩陣不同精度的觀測(cè)向量的值，并將其統(tǒng)計(jì)于表2。

表1 模擬真值統(tǒng)計(jì)

表2 僅含有偶然誤差模擬觀測(cè)值統(tǒng)計(jì)

對(duì)表2中的模擬觀測(cè)數(shù)據(jù)分別采用LS、WTLS、RTLS進(jìn)行計(jì)算，所得待估參數(shù)統(tǒng)計(jì)見表3。

表3 參數(shù)估值統(tǒng)計(jì)

由表3可以看出，當(dāng)系數(shù)矩陣含有誤差時(shí)，采用WTLS和RTLS參數(shù)估計(jì)精度高于LS。由于觀測(cè)值僅含有偶然誤差，在采用RTLS進(jìn)行迭代的過程中，權(quán)函數(shù)保持不變，因此，RTLS和WTLS所得結(jié)果相同。

3.2 觀測(cè)值含有偶然誤差和粗差

為了驗(yàn)證RTLS的穩(wěn)健性，在表2中第6組模擬觀測(cè)值中x和z分別添加2.0和1.0的粗差，其余組不變，然后采用LS、WTLS和RTLS對(duì)平面方程系數(shù)進(jìn)行估計(jì)，所得結(jié)果統(tǒng)計(jì)見表4。

由表4可知，當(dāng)觀測(cè)值混入粗差時(shí)，RTLS估計(jì)效果明顯好于LS估計(jì)和WTLS估計(jì)。這是由于在RTLS迭代過程中，含粗差觀測(cè)值的權(quán)函數(shù)元素的值會(huì)越來越小，并逐步趨近于0，不含粗差觀測(cè)值的權(quán)函數(shù)元素變化不大。因此，可以通過觀察迭代過程中權(quán)函數(shù)元素?cái)?shù)值的變化來確定粗差所在位置。在本次試驗(yàn)中，僅x2、x6、y3、y4和z6的權(quán)函數(shù)值發(fā)生變化，圖1為迭代過程中權(quán)函數(shù)值變化曲線圖。

表4 參數(shù)估值統(tǒng)計(jì)

圖1 權(quán)函數(shù)值變化曲線圖

由圖1可知，含有粗差的x6和z6觀測(cè)值權(quán)函數(shù)數(shù)值變化較快，隨著迭代過程，其值越來越小。平差過程就是誤差分配的過程。本方法采用平差后的殘差來確定權(quán)函數(shù)，當(dāng)觀測(cè)值混入粗差時(shí)必然會(huì)影響到其他不含粗差的觀測(cè)值，因此，也就解釋了為何不含粗差的x2、y3和y4觀測(cè)的權(quán)函數(shù)也會(huì)隨著迭代過程發(fā)生變化。

4 結(jié)論與不足

EIV模型同時(shí)考慮觀測(cè)向量和系數(shù)矩陣的誤差，因此自提出以來便得到各個(gè)領(lǐng)域廣泛的研究。針對(duì)EIV模型的估計(jì)方法統(tǒng)稱為整體最小二乘解。整體最小二乘經(jīng)過幾十年的發(fā)展，研究者們已經(jīng)提出了各種各樣的解法。然而這些解法主要都是考慮EIV模型中的偶然誤差，關(guān)于EIV模型的可靠性理論研究相對(duì)薄弱。本文針對(duì)EIV模型混入粗差的情況，導(dǎo)出穩(wěn)健整體最小二乘迭代解法(RTLS)。實(shí)驗(yàn)表明，RTLS通過觀察迭代過程中權(quán)函數(shù)元素?cái)?shù)值的變化來確定粗差所在位置，且其參數(shù)估值受粗差影響較小，具有穩(wěn)健性。然而，本文導(dǎo)出的RTLS解法是一種驗(yàn)后方法，其中權(quán)函數(shù)是根據(jù)每次迭代平差后的殘差來確定的，那些不含粗差的觀測(cè)值在平差過程中必然會(huì)受到含粗差的觀測(cè)值的影響。因此，如何在驗(yàn)前對(duì)平差模型進(jìn)行可靠性檢驗(yàn)還有待進(jìn)一步的研究與探索。

[1]陶本藻，邱衛(wèi)寧.誤差理論與測(cè)量平差[M].武漢：武漢大學(xué)出版社，2012.

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

[3]邱衛(wèi)寧，陶本藻，姚宜斌，等.測(cè)量數(shù)據(jù)處理理論與方法[M].武漢: 武漢大學(xué)出版社，2008.

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

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

[6]NEITZEL F.Generalization of Total Least-Squares on Example of Unweighted and Weighted 2D Similarity Transformation [J].Journal of Geodesy, 2010, 84(12):751-762.

[7]李德仁，袁修孝.誤差處理與可靠性理論[M].武漢：武漢大學(xué)出版社，2002:243-255.

[8]CHOI Y J, KIM J Y, SUMG K M.A Robust Algorithm of Total Least Squares Method[J].IEICETRANS,Fundamentals,1997,7(E80-A):1336-1339.

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

[10]VAHID MAHBOUB.On weighted total least-squares for geodetic transformations[J].Journal of Geodesy, 2012,86(5): 359-367.

[11]CHI LUN CHENG, MASTRONDRDL N, PALGE C.Total least squares and errors-in-variables moseling [J].Computational Statistics and Data Analysis,2007,52:1076-1079.

[12]V MAHBOUB, A R AMIRI-SIMKOOEI, M A SHARIFI.Iteratively reweighted total least squares: a robust estimation in errors-in-variables models[J].Survey Review,2013,45(329):92-99.

[13]KRARUP T, JUHL J, KUBIK K.Gotterdammerung over least squares adjustment[J].Proceedings of 14th Congress of the International Society of Photogrammetry,1980,B3: 369-378.

[14]卜長(zhǎng)江，羅躍生.矩陣論[M].哈爾濱：哈爾濱工程大學(xué)出版社，2007.

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

Robust estimation based on EIV model

CHU Bin, FAN Dong-ming, LIU Bo ,QIN Ning

(School of Earth Sciences and Environmental Engineering,Southwest Jiaotong University, Chengdu 610031，China)

EIV (error-in-variables) model has been widely used sin ce it is proposed as taking both the error of the observation vector and the coefficient matrix into account.However, the least squares solution for the EIV model which is called total least squares(TLS) assuming observations only contain accidental error, when there is a gross error in the observations,the solution is not optimal.By selecting an appropriate weight function, combined with the weighted total least squares (WTLS), a robust estimate called robust total least squares (RTLS) is proposed based on the EIV model.The linear fitting experiments show that the proposed method can locate gross errors, and the estimated amount is less affected by gross error, with robustness.

EIV model; total least squares(TLS); robust estimate; robust total least squares(RTLS); liner fitting

2013-08-11

楚 彬(1990-)，男，碩士研究生.

P207

：A

：1006-7949(2014)09-0017-04