線(xiàn)性EIV模型的TLS估計(jì)及其典型應(yīng)用

周擁軍，鄧才華

(1. 中南大學(xué) 有色金屬成礦預(yù)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室 長(zhǎng)沙 410083；2. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083)

現(xiàn)代測(cè)量平差與數(shù)據(jù)處理理論是以高斯?馬爾柯夫(Gauss-Markov, GM)模型為核心，通過(guò)該模型在不同層面上的拓展形成的若干新理論、新方法[1?2]。GM模型僅考慮觀測(cè)值的誤差，而認(rèn)為系數(shù)矩陣沒(méi)有誤差，而在實(shí)際應(yīng)用中，大多數(shù)情況下系數(shù)矩陣是有誤差的。EIV模型是指觀測(cè)值和參數(shù)均包含誤差的模型，針對(duì)線(xiàn)性EIV模型，GOLUB等[3]最早提出了總體最小二乘估計(jì)方法，MARKOVSKY 等[4?5]對(duì)擴(kuò)展 TLS問(wèn)題及其算法進(jìn)行了深入的研究，SCHAFFRIN 等[6?9]和NEITZEl[10]將TLS估計(jì)引入到測(cè)量數(shù)據(jù)處理中，國(guó)內(nèi)測(cè)繪領(lǐng)域的學(xué)者也從理論、應(yīng)用等方面對(duì)TLS進(jìn)行了深入研究[11?14]。但TLS用于解算線(xiàn)性GM模型和EIV模型估計(jì)方法的差異，各種擴(kuò)展TLS模型解算具體問(wèn)題的差別等尚未有深入的研究。

為方便總體最小二乘方法在測(cè)量數(shù)據(jù)處理中的應(yīng)用，這里將系數(shù)矩陣的誤差分為兩類(lèi)，一類(lèi)是線(xiàn)性化造成的設(shè)計(jì)矩陣元素的誤差，經(jīng)典GM模型大多是非線(xiàn)性的觀測(cè)方程在取參數(shù)近似值后得到的線(xiàn)性方程，參數(shù)初值對(duì)系數(shù)矩陣的取值有直接影響，這時(shí)系數(shù)矩陣包含的不僅僅是誤差而是包含模型誤差，這類(lèi)問(wèn)題在大多數(shù)需要線(xiàn)性化的測(cè)量平差問(wèn)題中大量存在，而且需要迭代計(jì)算。另一類(lèi)是設(shè)計(jì)矩陣元素是只包含觀測(cè)值，如線(xiàn)性回歸、直接線(xiàn)性變換、曲線(xiàn)擬合等，很顯然設(shè)計(jì)矩陣含有誤差。本文作者比較了GM模型與線(xiàn)性EIV模型和估計(jì)方法的區(qū)別，并給出了兩類(lèi)問(wèn)題的算例，旨在推廣TLS估計(jì)在測(cè)量數(shù)據(jù)處理中的應(yīng)用。

1 線(xiàn)性EIV模型與GM模型的比較

測(cè)量平差領(lǐng)域常用的線(xiàn)性GM模型可表示為

式中： L ∈ Rm×1表示觀測(cè)值向量， Δ ∈ Rm×1表示觀測(cè)值向量的真誤差， X ∈Rn×1表示待估參數(shù)，A∈Rm×n表示系數(shù)矩陣，其中 m、n分別表示觀測(cè)值和未知參數(shù)的個(gè)數(shù)，且滿(mǎn)足m≥n。 σ02表示單位權(quán)方差，P表示觀測(cè)值的權(quán)矩陣。經(jīng)典測(cè)量平差問(wèn)題是以最小二乘準(zhǔn)則ΔTPΔ=min為代價(jià)函數(shù)的參數(shù)估計(jì)，可以證明最小二乘準(zhǔn)則和最大似然準(zhǔn)則是一致的，具有無(wú)偏，方差最小等統(tǒng)計(jì)特征。

線(xiàn)性EIV(Errors-in-variables)模型除了考慮觀測(cè)值的誤差外，還考慮矩陣A的誤差，其數(shù)學(xué)模型可表示為

式中：EA表示觀測(cè)矩陣 A的誤差向量，A~表示系數(shù)矩陣的真值， eA= vec(EA)∈ Rmn×1表示矩陣EA按列向量化后得到的向量，PA∈Rmn×mn表示 eA的權(quán)陣。對(duì)于線(xiàn)性 EIV模型應(yīng)采用總體最小二乘估計(jì)，文獻(xiàn)[3]中已證明總體最小二乘估計(jì)是EIV模型的最大似然估計(jì)，其目標(biāo)函數(shù)可表示為

為簡(jiǎn)化表達(dá)，先假設(shè)所有誤差服從獨(dú)立、等方差的正態(tài)分布，則權(quán)矩陣為單位矩陣，式(3)可表示為

引入 Eckart-Young-Mirsky矩陣近似定理：矩陣A∈Rm×n，其秩 rank(A)=r的 svd分解可表示為 A=其中，σi為矩陣A的奇異值，U、V為正交矩陣，其向量形式為 U =(u1u2…um)，V=(v1v2…vn)，uj和vj分別為矩陣A的奇異值σj對(duì)應(yīng)的左、右向量。若k≤r，，并存在以下關(guān)系：

若不考慮參數(shù)之間的約束則rank(A)=n，即系數(shù)矩陣列滿(mǎn)秩，欲使0有非零解，則必有rank (A)＜n+1，根據(jù) Eckart-Young-Mirsky矩陣近似定理，取rank ( A?)=n作為原矩陣的近似，則有：

并得到待求參數(shù)為

以上推導(dǎo)表明，經(jīng)典 GM 模型可以擴(kuò)展如式(6)所示的線(xiàn)性EIV模型，其對(duì)應(yīng)的估計(jì)方法也由傳統(tǒng)的LS估計(jì)方法變?yōu)門(mén)LS估計(jì)，對(duì)應(yīng)的擴(kuò)展參數(shù)解X的解為經(jīng)奇異值分解后得到的最小特征值所對(duì)應(yīng)的右特征向量vn+1，由于=[X, ?1]，從而得到參數(shù)的估計(jì)值為

可以證明，它與常規(guī)的法方程求逆 X ? = ?(ATA?σI)?1ATL的估計(jì)結(jié)果是一致的，很顯然直接分解

n+1方法更簡(jiǎn)潔，得到參數(shù)的估計(jì)值后可進(jìn)一步得到的單位權(quán)方差：

2 TLS估計(jì)的擴(kuò)展

2.1 WTLS問(wèn)題

通過(guò)上面的分析可知，經(jīng)典GM模型可以方便地?cái)U(kuò)展為線(xiàn)性EIV模型，但目前諸多應(yīng)用并未考慮設(shè)計(jì)矩陣的權(quán)，所得到的總體二乘結(jié)果并非滿(mǎn)足統(tǒng)計(jì)上的最優(yōu)。欲確定權(quán)陣，首先要確定A的協(xié)方差陣C(A)，則C(A)∈Rmn×nm，若用A~表示A的真值，令A(yù)=A~+EA，根據(jù)協(xié)方差陣的定義：

C(A)=E[(EA)ij(EA)kl] (12)式中：i，j，k，l表示對(duì)應(yīng)的矩陣元素。若考慮協(xié)方差陣的一般形式，勢(shì)必增加計(jì)算的工作量，MARKOVSKY等[4?5]根據(jù)不同的應(yīng)用實(shí)例對(duì)方差陣進(jìn)行了簡(jiǎn)化，若矩陣A按行獨(dú)立，且每行具有相同的方差陣，則這類(lèi)問(wèn)題稱(chēng)為廣義總體最小二乘問(wèn)題(GTLS)。若矩陣A按行獨(dú)立，但各行的方差陣不同，則稱(chēng)為按行的加權(quán)總體最小二乘問(wèn)題(Row-wise weighted total least-squares problem，RWTLS)。根據(jù)權(quán)陣的結(jié)構(gòu)，可以將加權(quán)總體最小二乘問(wèn)題表示為表 1的形式，除廣義最小二乘可以采用類(lèi)似簡(jiǎn)單TLS問(wèn)題的直接分解算法外，其它的WTLS均需要迭代[4]。

測(cè)量平差中的GM模型，系數(shù)矩陣A是獨(dú)立的觀測(cè)方程經(jīng)線(xiàn)性化得到的，即：

式中：

由于觀測(cè)值間相互獨(dú)立，矩陣A的各行相互獨(dú)立，但每行的元素是相關(guān)的，因此屬于 Row-wise WTLS問(wèn)題，用αA表示A的第α行(],1[m∈α)，其對(duì)應(yīng)的方差陣如下：

C(X)表示未知參數(shù)的方差陣由下式確定：

從而得到擴(kuò)展參數(shù)的方差陣為C ( X)，進(jìn)而簡(jiǎn)單總體最小二乘問(wèn)題變換為加權(quán)最小二乘問(wèn)題：于權(quán)陣與參數(shù)相關(guān)，需要迭代求解，具體的計(jì)算步驟如下[4]：

1) 根據(jù)給定的參數(shù)初值列出誤差方程；

2) 給定權(quán)陣PA=I按簡(jiǎn)單 TLS方法得到參數(shù)初值；

3) 根據(jù)誤差傳播率得到參數(shù)的方差陣 C(X)，并逐行得到各元素的方差陣C(Aα)；

4) 根據(jù)加權(quán)則解為將S進(jìn)行SVD分解后最小特征值所對(duì)應(yīng)的特征向量；

5) 將得到的參數(shù)與初值比較，如果小于某一閾值，則計(jì)算完成，否則回到式(3)進(jìn)行迭代計(jì)算。

TLS可以通過(guò)奇異值分解算法實(shí)現(xiàn)，而WTLS由

2.2 混合OLS-TLS問(wèn)題

測(cè)量中的權(quán)矩陣為方差陣的逆矩陣，而系數(shù)矩陣A的某些元素可能為常數(shù)，因此這些項(xiàng)對(duì)應(yīng)的方差分量為零，從而導(dǎo)致無(wú)法通過(guò)方差陣求逆得到權(quán)陣，此時(shí)需要采用混合最小二乘方法求解，這里先不考慮A的權(quán)，將式(1)變?yōu)?/p>

表1 權(quán)陣的結(jié)構(gòu)及對(duì)應(yīng)的總體最小二乘問(wèn)題Table 1 Weighted matrix structure and its corresponding TLS problems

式中：系數(shù)矩陣 A1無(wú)誤差，而 A2含有誤差，則A=[A , A , L]，X = [XT, XT, ?1]T，X =[X2, ?1]，將A

1 21 22

進(jìn)行QR分解，得到：

可以根據(jù)R22X2=0解出 X2，然后回代到方程R11X1+R12X2=0中，解出 X1的值，然后根據(jù)前面的方法進(jìn)行精度評(píng)定。

2.3 CTLS問(wèn)題

這里僅考慮線(xiàn)性約束條件，至于非線(xiàn)性約束條件可以根據(jù)參數(shù)初值線(xiàn)性化后得到線(xiàn)性約束條件，線(xiàn)性條件一般可表達(dá)為

CTLS問(wèn)題可以采用迭代解法[4]或直接分解算法[15]，具體的解算過(guò)程本文不再列出。

3 算例

選擇兩個(gè)典型算例，一個(gè)是典型的大地測(cè)量控制網(wǎng)平差，需要給定參數(shù)初值轉(zhuǎn)化為誤差方程，然后根據(jù)間接平差法求解(以下簡(jiǎn)稱(chēng)LS估計(jì))，然后將誤差方程拓展為線(xiàn)性EIV模型，用TLS方法求解，并比較二者在參數(shù)初值、迭代次數(shù)、估計(jì)精度等的差別；二是曲線(xiàn)擬合的算例，其設(shè)計(jì)矩陣有誤差，但不受參數(shù)初值的影響，比較分析各種估計(jì)方法。

3.1 測(cè)量控制網(wǎng)平差算例

圖1所示為一測(cè)角網(wǎng)，A、B、C 3點(diǎn)的坐標(biāo)分別為(8 864.530，5 392.580)、(13 615.220，10 189.470)、(6 520.120、14 399.300)，觀測(cè)數(shù)據(jù)見(jiàn)表1。為了比較初值對(duì)估計(jì)結(jié)果的影響，選擇了兩組初值，分別采用經(jīng)典最小二乘方法、總體最小二乘方法和加權(quán)總體最小二乘方法進(jìn)行迭代計(jì)算，?。?0?6作為迭代終止條件，計(jì)算的結(jié)果見(jiàn)表 3。計(jì)算結(jié)果表明：當(dāng)初值接近估計(jì)值時(shí)，3種方法的迭代次數(shù)和最終估值基本相等，而經(jīng)典大地測(cè)量網(wǎng)的初值往往與估值較接近，說(shuō)明迭代LS估計(jì)和TLS估計(jì)在計(jì)算速度和精度上并無(wú)明顯改進(jìn)。若初值偏離估值較遠(yuǎn)，迭代TLS估計(jì)要略快于 LS估計(jì)的收斂速度，在估計(jì)結(jié)果上也僅有微小差異，說(shuō)明初值和估計(jì)方法對(duì)結(jié)果影響不大。但需要指出的是，簡(jiǎn)單TLS在設(shè)計(jì)矩陣的所有元素都滿(mǎn)足獨(dú)立等精度分布的前提下才是最優(yōu)估計(jì)[4]，而這一隨機(jī)假設(shè)在實(shí)際問(wèn)題中基本不成立，導(dǎo)致簡(jiǎn)單TLS估計(jì)多為有偏估計(jì)。當(dāng)函數(shù)模型較復(fù)雜、設(shè)計(jì)矩陣病態(tài)時(shí)，TLS估計(jì)和LS估計(jì)的結(jié)果就會(huì)有較大差異，故簡(jiǎn)單TLS僅適用于近似平差計(jì)算。在嚴(yán)密平差時(shí)，需采用WTLS方法或經(jīng)典平差方法，但WTLS方法需要根據(jù)觀測(cè)值的誤差傳播率得到設(shè)計(jì)矩陣元素的方差，以設(shè)計(jì)矩陣元素的最小二乘準(zhǔn)則為目標(biāo)函數(shù)進(jìn)行平差計(jì)算，這導(dǎo)致計(jì)算較間接平差方法復(fù)雜，因此，經(jīng)典控制網(wǎng)平差從計(jì)算效率和精度而言并不適宜采用WTLS方法。

圖1 某大地測(cè)量控制網(wǎng)Fig. 1 A typical geodetic control network

3.2 曲線(xiàn)擬合算例

便于和為真值比較，這里采用模擬算例，設(shè)二次曲線(xiàn)的方程為y=ax2+bx+c，參數(shù)真值分別為：3、5、10，取x為 [0,10]區(qū)間內(nèi)產(chǎn)生的10個(gè)隨機(jī)數(shù)，并按以上方程得到y(tǒng)的值，然后在x和y上均疊加方差為0.01的正態(tài)分布誤差，分別用LS、TLS、OLS-TLS、WTLS估計(jì)參數(shù)，做10 000次模擬實(shí)驗(yàn)，然后比較各種估計(jì)方法的均方誤差得到的計(jì)算結(jié)果見(jiàn)表4，表明TLS和OLS-TLS的計(jì)算效果甚至比LS還差，主要原因是TLS雖然考慮了設(shè)計(jì)矩陣的誤差，但假設(shè)矩陣各元素的誤差是獨(dú)立等精度分布的，而本實(shí)例中各元素的誤差明顯不相等，因此造成了TLS估計(jì)效果不理想，而WTLS由于考慮了各元素的統(tǒng)計(jì)特性，因而估計(jì)效果最佳。作者還采用變換各參數(shù)數(shù)量級(jí)的方法來(lái)模擬，發(fā)現(xiàn)簡(jiǎn)單TLS和LS方法各有優(yōu)劣，但WTLS的估計(jì)效果仍為最佳，由于篇幅原因在此不再列出。

表2 觀測(cè)數(shù)據(jù)Table 2 Observed data

表3 平差結(jié)果Table 3 Results estimated by LS, TLS and WTLS

表4 LS、TLS、OLS-TLS和 WTLS方法的參數(shù)均方誤差比較Table 4 Roots of mean square estimated by LS, TLS,OLS-TLS and WTLS

4 結(jié)論和建議

1) 線(xiàn)性GM模型能方便地轉(zhuǎn)換為線(xiàn)性EIV模型并用TLS方法求解。

2) 對(duì)于線(xiàn)性EIV模型，簡(jiǎn)單TLS估計(jì)僅在設(shè)計(jì)矩陣元素服從獨(dú)立等方差分布時(shí)才是最優(yōu)估計(jì)，因此，簡(jiǎn)單TLS方法的估計(jì)精度不一定優(yōu)于LS方法的估計(jì)精度，簡(jiǎn)單TLS估計(jì)一般適用于近似平差計(jì)算。

3) 對(duì)于不等精度觀測(cè)值問(wèn)題，若要獲取高精度的平差結(jié)果，宜采用WTLS方法。雖然經(jīng)典控制網(wǎng)平差的函數(shù)模型可經(jīng)線(xiàn)性化得到線(xiàn)性GM模型，而后轉(zhuǎn)換為EIV模型并采用WTLS平差方法得到與經(jīng)典方法精度相當(dāng)?shù)钠讲罱Y(jié)果，但計(jì)算過(guò)程較間接平差復(fù)雜，因此建議仍采用經(jīng)典間接平差方法。而對(duì)于曲線(xiàn)擬合等以坐標(biāo)為觀測(cè)值的平差問(wèn)題，采用WTLS估計(jì)相對(duì)簡(jiǎn)便，本文作者將做進(jìn)一步的研究。

[1] 朱建軍, 宋迎春. 現(xiàn)代測(cè)量平差與數(shù)據(jù)處理理論的進(jìn)展[J].工程勘察, 2009, 37(12): 1?5.

ZHU Jian-jun, SONG Ying-chun. Progress of modern surveying adjustment and theory of data processing [J]. Geotechnical Investigation & Surveying, 2009, 37(12): 1?5.

[2] 歐陽(yáng)文森, 朱建軍. 經(jīng)典平差模型的擴(kuò)展[J]. 測(cè)繪學(xué)報(bào), 2009,38(1): 12?15.

OUYANG Wen-sen, ZHU Jian-jun. Expanding of classical surveying adjustment model [J]. Acta Geodaeticaet Cartographic Sinica, 2009, 38(1): 12?15.

[3] GOLUB G H, VAN LOAN C F. An analysis of the total least squares problem [J]. SIAM J Numer Anal, 1980, 17: 883?893.

[4] MARKOVSKY I, VAN HUFFEL S. Overview of total least squares methods [J]. Signal Process, 2007, 87(10): 2283?2302.

[5] MARKOVSKY I, RASTELLO M L, PREMOLI P, KUKUSIT A,van HUFFEL S. The element-wise weighted total least squares problem [J]. Computational Statistics & Data Analysis, 2006,50(1): 181?209.

[6] SCHAFFRIN B, FELUS Y A. On total least-squares adjustment with constraints [C]// Windows on the Future of Geodesy,IAG-Symp. Springer, Berlin, 2005, 128: 417?421.

[7] SCHAFFRIN B. A note on constrained total least squares estimation [J]. Linear Algebra Application, 2006, 417(1): 245?258.

[8] SCHAFFRIN B, WIESER A. On weighted total least-squares adjustment for linear regression [J]. Journal of Geodesy, 2008,82(6): 373?383.

[9] 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.

[10] NEITZEL F. Generalization of total least-squares on example of unweighted and weighted 2D similarity transformation [J].Journal of Geodesy, 2010, 84(12): 373?383.

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

LU Tie-ding, TAO Ben-zao, ZHOU Shi-jian. Modeling and algorithm of linear regression based on total least squares [J].Geomatics and Information Science of Wuhan University, 2008,33(5): 504?507.

[12] 陸 玨, 陳 義, 鄭 波. 總體最小二乘方法在三維坐標(biāo)轉(zhuǎn)換中的應(yīng)用[J]. 大地測(cè)量與地球動(dòng)力學(xué), 2008, 28(5): 77?81.

LU Jue, CHEN Yi, ZHENG Bo. Applying total least squares to three dimensional datum transformation [J]. Journal of Geodesy and Geodynamics, 2008, 28(5): 77?81.

[13] 袁振超, 沈云中, 周澤波. 病態(tài)總體最小二乘模型的正則化算法[J]. 大地測(cè)量與地球動(dòng)力學(xué), 2009, 29(2): 131?134.

YUAN Zhen-chao, SHEN Yun-zhong, ZHOU Ze-bo.Regularized total least-squares solution to ill-posed error-invariable model [J]. Journal of Geodesy and Geodynamics, 2009,29(2): 131?134.

[14] 孔 建, 姚宜斌, 吳 寒. 整體最小二乘的迭代解法[J]. 武漢大學(xué)學(xué)報(bào): 信息科學(xué)版, 2010, 35(6): 711?714.

KONG Jian, YAO Yi-bin, WU Han. Iterative method for total least-squares [J]. Geomatics and Information Science of Wuhan University, 2010, 35(6): 711?714.

[15] DOWLING E M, DEGROAT R D, LINEBARGER D A. Total least squares with linear constraints [C]// IEEE, International Conference on Acoustics, Speech, and Signal Processing, 1992:341?347.