線性EIV模型的TLS估計及其典型應用

周擁軍，鄧才華

(1. 中南大學 有色金屬成礦預測教育部重點實驗室 長沙 410083；2. 中南大學 地球科學與信息物理學院，長沙 410083)

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

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

1 線性EIV模型與GM模型的比較

測量平差領域常用的線性GM模型可表示為

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

線性EIV(Errors-in-variables)模型除了考慮觀測值的誤差外，還考慮矩陣A的誤差，其數學模型可表示為

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

為簡化表達，先假設所有誤差服從獨立、等方差的正態分布，則權矩陣為單位矩陣，式(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對應的左、右向量。若k≤r，，并存在以下關系：

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

并得到待求參數為

以上推導表明，經典 GM 模型可以擴展如式(6)所示的線性EIV模型，其對應的估計方法也由傳統的LS估計方法變為TLS估計，對應的擴展參數解X的解為經奇異值分解后得到的最小特征值所對應的右特征向量vn+1，由于=[X, ?1]，從而得到參數的估計值為

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

n+1方法更簡潔，得到參數的估計值后可進一步得到的單位權方差：

2 TLS估計的擴展

2.1 WTLS問題

通過上面的分析可知，經典GM模型可以方便地擴展為線性EIV模型，但目前諸多應用并未考慮設計矩陣的權，所得到的總體二乘結果并非滿足統計上的最優。欲確定權陣，首先要確定A的協方差陣C(A)，則C(A)∈Rmn×nm，若用A~表示A的真值，令A=A~+EA，根據協方差陣的定義：

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

測量平差中的GM模型，系數矩陣A是獨立的觀測方程經線性化得到的，即：

式中：

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

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

從而得到擴展參數的方差陣為C ( X)，進而簡單總體最小二乘問題變換為加權最小二乘問題：于權陣與參數相關，需要迭代求解，具體的計算步驟如下[4]：

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

2) 給定權陣PA=I按簡單 TLS方法得到參數初值；

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

4) 根據加權則解為將S進行SVD分解后最小特征值所對應的特征向量；

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

TLS可以通過奇異值分解算法實現，而WTLS由

2.2 混合OLS-TLS問題

測量中的權矩陣為方差陣的逆矩陣，而系數矩陣A的某些元素可能為常數，因此這些項對應的方差分量為零，從而導致無法通過方差陣求逆得到權陣，此時需要采用混合最小二乘方法求解，這里先不考慮A的權，將式(1)變為

表1 權陣的結構及對應的總體最小二乘問題Table 1 Weighted matrix structure and its corresponding TLS problems

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

1 21 22

進行QR分解，得到：

可以根據R22X2=0解出 X2，然后回代到方程R11X1+R12X2=0中，解出 X1的值，然后根據前面的方法進行精度評定。

2.3 CTLS問題

這里僅考慮線性約束條件，至于非線性約束條件可以根據參數初值線性化后得到線性約束條件，線性條件一般可表達為

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

3 算例

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

3.1 測量控制網平差算例

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

圖1 某大地測量控制網Fig. 1 A typical geodetic control network

3.2 曲線擬合算例

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

表2 觀測數據Table 2 Observed data

表3 平差結果Table 3 Results estimated by LS, TLS and WTLS

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

4 結論和建議

1) 線性GM模型能方便地轉換為線性EIV模型并用TLS方法求解。

2) 對于線性EIV模型，簡單TLS估計僅在設計矩陣元素服從獨立等方差分布時才是最優估計，因此，簡單TLS方法的估計精度不一定優于LS方法的估計精度，簡單TLS估計一般適用于近似平差計算。

3) 對于不等精度觀測值問題，若要獲取高精度的平差結果，宜采用WTLS方法。雖然經典控制網平差的函數模型可經線性化得到線性GM模型，而后轉換為EIV模型并采用WTLS平差方法得到與經典方法精度相當的平差結果，但計算過程較間接平差復雜，因此建議仍采用經典間接平差方法。而對于曲線擬合等以坐標為觀測值的平差問題，采用WTLS估計相對簡便，本文作者將做進一步的研究。

[1] 朱建軍, 宋迎春. 現代測量平差與數據處理理論的進展[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] 歐陽文森, 朱建軍. 經典平差模型的擴展[J]. 測繪學報, 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] 魯鐵定, 陶本藻, 周世健. 基于整體最小二乘法的線性回歸建模和解法[J]. 武漢大學學報: 信息科學版, 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] 陸 玨, 陳 義, 鄭 波. 總體最小二乘方法在三維坐標轉換中的應用[J]. 大地測量與地球動力學, 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] 袁振超, 沈云中, 周澤波. 病態總體最小二乘模型的正則化算法[J]. 大地測量與地球動力學, 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]. 武漢大學學報: 信息科學版, 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.