方差分量估計的抗差總體最小二乘算法

陶葉青，蔡安寧，楊 娟

(淮陰師范學院 城市與環(huán)境學院，江蘇 淮安 223300)

為解決變量中誤差的模型參數(shù) (errors-in-variables, EIV)估計問題，最小二乘約束準則(least squares estimation, LS)被拓展為總體最小二乘約束準則(total least squares estimation, TLS)[1]。應用總體最小二乘約束準則建立模型參數(shù)估計算法受到研究人員的關(guān)注，此算法可以使總體最小二乘參數(shù)估計模型被應用到坐標相似變換[2-9]、高程異常擬合[10]、大地測量反演[11]等測量數(shù)據(jù)處理中。

研究人員在討論模型參數(shù)估計算法的同時，將經(jīng)典平差模型中存在的粗差處理問題引入到總體最小二乘模型中進行討論[12-18]。文獻[12]應用粗差識別理論闡述總體最小二乘模型粗差處理的方法，但是算法沒有通過實例進行驗證。基于粗差調(diào)節(jié)理論，應用等價權(quán)函數(shù)建立的抗差總體最小二乘算法(robust total least squares, RTLS)是處理含有粗差的EIV模型的主要方法：文獻[13]基于Huber權(quán)函數(shù)，應用高斯-赫爾墨特模型(Gauss-Helmert model, G-H)建立總體最小二乘抗差迭代算法，并將其應用至空間坐標系統(tǒng)轉(zhuǎn)換；隨后，文獻[14]～[18]分別基于不同的權(quán)函數(shù)，應用極值函數(shù)建立含有粗差的總體最小二乘模型抗差估計算法。

研究人員采用經(jīng)典平差模型的粗差處理思路解決總體最小二乘模型的粗差處理問題，需要注意的是含有粗差的總體最小二乘模型具有顯著的特征，其研究思路不能完全趨同于最小二乘模型。此類問題引起研究人員的關(guān)注：總體最小二乘模型系數(shù)矩陣與觀測向量可能由不同類型的觀測元素組成，直接通過觀測值的改正數(shù)調(diào)節(jié)權(quán)值的抗差算法，在文獻[19]～[21]得到改進。文獻[19]借助標準化殘差，文獻[20]借助敏感度分析構(gòu)造統(tǒng)計量，增加抗差總體最小二乘算法的通用性。同時，總體最小二乘模型是非線性模型，在有限觀測元素條件下，得到的單位權(quán)中誤差是有偏的[3]。文獻[7]建立了單位權(quán)中誤差無偏估計的迭代算法，但是沒有顧及含有粗差的觀測值對估值的影響。此外，當模型的系數(shù)矩陣與觀測向量由不同類型的觀測元素組成時，由相同的中誤差確定權(quán)函數(shù)的域值是不合理的。為合理確定等價權(quán)函數(shù)中分段函數(shù)的域值，中位數(shù)法被應用在抗差總體最小二乘估計[21]?，F(xiàn)有研究成果是假定隨機模型是精確、已知的，缺乏當隨機模型存在偏差對抗差算法影響的探討。

為顧及隨機模型誤差對總體最小二乘抗差估計的影響，有效應用等價權(quán)函數(shù)實現(xiàn)抗差估計，本文討論采用方差分量估計(variance component estimation, VCE)[22]修正隨機模型，通過最小二乘方差分量估計(least squares variance component estimation, LS-VCE)的總體最小二乘參數(shù)估計迭代算法[23]，分別計算系數(shù)矩陣與觀測向量的中誤差。在應用權(quán)函數(shù)對觀測元素進行迭代定權(quán)的同時，對隨機模型進行修證，分類確定不同類型觀測值對應權(quán)函數(shù)的域值。

1 抗差總體最小二乘估計

總體最小二乘準則能夠同時對系數(shù)矩陣與觀測向量中含有的誤差附加約束：

(1)

式中：Δ為觀測向量L中含有的誤差向量，EB為系數(shù)矩陣B中含有的誤差矩陣，x為參數(shù)向量，b=vec(EB)(vec(·)為矩陣拉直符號，表示將誤差矩陣按列拉直，得到列向量)；P，PB為觀測向量與系數(shù)矩陣的權(quán)矩陣，其中PB可以由協(xié)因數(shù)陣的Kron積得到，具體定義過程請參見文獻[3]。需要指出的是，對于partial-EIV模型[24]系數(shù)矩陣中常數(shù)元素的定權(quán)問題，可以通過協(xié)因數(shù)陣的定義來解決。

(2)

(3)

式中：A0為拉格朗日函數(shù)關(guān)于參數(shù)向量的偏導數(shù)，x0為參數(shù)向量的初值，ω0=L-Bx0，B10為拉格朗日函數(shù)關(guān)于誤差向量Δ的偏導數(shù)，B20為關(guān)于誤差向量b的偏導數(shù)，N是關(guān)于B10與B20的函數(shù)，λ為拉格朗日乘數(shù)，Q1為參數(shù)向量的協(xié)因數(shù)陣，Q2為系數(shù)矩陣的協(xié)因數(shù)陣。模型參數(shù)的具體求解方法請參見文獻[15]。根據(jù)系數(shù)矩陣與觀測向量中殘差的估值，應用等價權(quán)函數(shù)對觀測元素進行重新定權(quán)，以IGGⅡ權(quán)函數(shù)為例[25]

(4)

抗差總體最小二乘參數(shù)估計的研究思路與經(jīng)典平差模型的抗差估計研究思路趨于一致，其優(yōu)點是理論基礎(chǔ)成熟、算法容易實現(xiàn)，但是忽略了總體最小二乘模型與最小二乘模型之間的差異。在總體最小二乘模型中，當系數(shù)矩陣與觀測向量中的元素不屬于同一類觀測值時，兩者在精度上往往不屬于同一量級，并且隨機模型的定義存在偏差。此時，通過單一的中誤差統(tǒng)一確定觀測值的權(quán)，這樣的方法顯然是不合理的。有學者提出“附有相對權(quán)比”[26]、“標準化殘差”[19]等方法來克服傳統(tǒng)算法存在的問題，但是并不能夠消除隨機模型誤差對參數(shù)估計的影響。由此，論文探討應用方差分量估計建立總體最小二乘抗差估計算法。

2 基于方差分量估計的抗差算法

2.1 方差分量估計的最小二乘算法

方差分量估計是進行隨機模型驗后估計的有效方法，一直受到研究人員的關(guān)注，主要算法有赫爾默特方差估計、最小范數(shù)二次無偏估計(MINQUE)、最優(yōu)不變二次無偏估計(BIQUE)。最小二乘方差分量估計被證明優(yōu)于其它方差估計算法[22]，應用LS-VEC改正總體最小二乘隨機模型，在此基礎(chǔ)上建立抗差總體最小二乘參數(shù)估計的迭代算法。當不顧及系數(shù)矩陣中的誤差，EIV模型退化為高斯-馬爾可夫模型(Gauss-Markov model, G-M)

Δ=Bx-L.

(5)

設(shè)觀測向量L對應的協(xié)方差陣為D(l)，將其表示為觀測元素中誤差的函數(shù)

(6)

(7)

(8)

(9)

2.2 抗差總體最小二乘的迭代算法

EIV模型系數(shù)矩陣B對應的協(xié)方差陣為D(B)，將其表示為系數(shù)矩陣中觀測元素的中誤差函數(shù)

(10)

應用LS-VEC建立抗差總體最小二乘迭代算法：

2)應用式(6)、式(10)計算觀測向量與系數(shù)矩陣的協(xié)因數(shù)陣Ql，QB；

3)根據(jù)式(2)、式(3)計算EIV模型的總體最小二乘解，以及系數(shù)矩陣與觀測向量中誤差向量的估值，并用其對觀測向量與系數(shù)矩陣進行改正，由改正后的觀測向量、系數(shù)矩陣與調(diào)節(jié)后的協(xié)因數(shù)陣計算正交投影矩陣R，元素nij，li，以確定N1，l1，并且應用式(7)計算待估向量σk；

4)應用權(quán)函數(shù)對觀測元素進行分類定權(quán)。需要指出的是，當EIV模型的觀測向量與系數(shù)矩陣中的觀測元素屬于不同的精度類型，權(quán)函數(shù)的域值也可以采用標準化殘差[19]、敏感度分析[20]、中位數(shù)法[21]等進行分類確定；

5)迭代步驟(2)、(3)、(4)，直到參數(shù)估值收斂。

3 算例與分析

應用直線擬合對建立的算法進行驗證，當有n個觀測值時，將直線方程表示為

(11)

式中：[kb]T為模型待估參數(shù)向量。應用文獻[27]的數(shù)據(jù)作為實驗數(shù)據(jù)，見表1。

圖1 參數(shù)最優(yōu)估值的擬合精度

應用上述觀測數(shù)據(jù)對論文建立的算法進行驗證，分別在1號點的x坐標、3號點的y坐標中加入1個單位的誤差，模擬觀測值中含有的粗差；并且將求解參數(shù)的元素分為三組：一組(1，2，4，…，10)，二

組(2，3，4，…，10)，三組(1，2，3，…，10)，分別模擬模型系數(shù)矩陣中的元素含有粗差，模型觀測向量中的元素含有粗差，模型觀測向量與系數(shù)矩陣中的元素同時含有粗差的三種情況。應用基于等價權(quán)函數(shù)的抗差總體最小二乘傳統(tǒng)算法，不同觀測樣本求解的參數(shù)估值見表2。需要指出的是，三組樣本獲得的殘差改正數(shù)均小于權(quán)函數(shù)中第一個分段函數(shù)的域值，這顯然是和實際情況相違背的，權(quán)函數(shù)沒有對隨機模型進行調(diào)節(jié)。三組樣本擬合的精度見表3。

表2 模型參數(shù)的估值

基于方差分量估計的抗差算法，對系數(shù)矩陣元素與觀測向量元素分類確定權(quán)函數(shù)的域值，求解模型參數(shù)估值。經(jīng)過四次迭代計算，獲得參數(shù)的估值見表2，其擬合精度見表3，等價權(quán)函數(shù)能夠?qū)﹄S機模型進行有效的調(diào)節(jié)。

兩種算法擬合精度的誤差分布如圖2所示。

表3 不同算法擬合的模型參數(shù)精度

圖2 兩種算法擬合的誤差分布

實驗結(jié)果表明，應用方差分量估計建立的總體最小二乘抗差估計算法是可行的，建立算法的擬合精度優(yōu)于傳統(tǒng)總體最小二乘抗差算法。在進行的另一組模擬實驗中，系數(shù)矩陣與觀測向量中的元素精度存在較大差異，模擬的隨機模型存在偏差，獲得模型參數(shù)的精度更加顯著。

4 結(jié) 論

通過對現(xiàn)有總體最小二乘抗差估計研究成果的分析，歸納現(xiàn)有算法中存在的不足?？傮w最小二乘抗差估計的研究思路趨同于傳統(tǒng)算法，無法顧及EIV模型存在的系數(shù)矩陣與觀測向量中的元素不屬于同一精度量級引起的抗差估計問題，同時沒有顧及隨機模型誤差對基于權(quán)函數(shù)的TLS抗差估計算法影響的問題。應用最小二乘方差分量估計建立總體最小二乘抗差估計算法，修正隨機模型誤差對抗差估計的影響，建立的迭代算法能夠?qū)崿F(xiàn)對系數(shù)矩陣與觀測向量中的元素進行分類定權(quán)，提高應用權(quán)函數(shù)進行抗差估計的有效性。

[1] PEARSON K. On Lines and Planes of Closest Fit to Systems of Points in Space[J]. Phil Mag, 1901,2:559-572.

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

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

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

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

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

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

[8] 方興，曾文憲，劉經(jīng)南，等. 三維坐標轉(zhuǎn)換的通用整體最小二乘算法[J].測繪學報，2014， 43(11)：1139-1143.

[9] XING Fang. Weighted total least-squares with constraints: a universal formula for geodetic symmetrical transformations[J]. Journal of Geodesy, 2015, 89(5): 459-469.

[10] TAO Y Q, GAO J X, YAO Y F. TLS algorithm for GPS height fitting based on robust estimation[J]. Survey Review, 2014,46(336):184-188.

[11] 王樂洋，許才軍，魯鐵定.邊長變化反演應變參數(shù)的總體最小二乘方法[J].武漢大學學報(信息科學報)，2010，35(2)：181-184．

[12] SCHAFFRIN B, UZUN S. Errors-in-variables for mobile mapping algorithms in the presence of outliers[J]. Archives of Photogrammetry, 2011, 22: 377-387.

[13] 陳義，陸玨. 以三維坐標轉(zhuǎn)換為例解算穩(wěn)健總體最小二乘方法[J].測繪學報，2012，41(5)：715-722．

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

[15] 龔循強，李志林. 穩(wěn)健加權(quán)總體最小二乘法[J]. 測繪學報，2014，43(9)：888-894.

[16] 龔循強，李志林.一種利用IGGⅡ方案的穩(wěn)健混合總體最小二乘方法[J].武漢大學學報(信息科學報)，2014，39(4)：462-466．

[17] LU J, CHEN Y, LI B F,et al. Robust total least squares with reweighting iteration for three-dimensional similarity transformation[J]. Survey Review, 2014, 46(334): 28-36.

[18] 王彬，李建成，高井祥，等.抗差加權(quán)整體最小二乘模型的牛頓-高斯算法[J]. 測繪學報，2015，44(6)：602-608.

[19] WANG Bin, LI Jiancheng, LIU Chao. A robust weighted total least squares algorithm and its geodetic applications[J]. Stud. Geophys. Geod., 2016, 60: 177-194.

[20] 趙俊，歸慶明. 部分變量誤差模型的整體抗差最小二乘估計[J]. 測繪學報，2016，45(5)：552-559.

[21] 陶葉青，高井祥，姚一飛．基于中位數(shù)法的抗差總體最小二乘估計[J]．測繪學報，2015，45(3)：297-301．

[22] TEUNISSEN P J G, AMIRI-SIMKOOEI A R. Least-squares variance component estimation[J]. Journal of Geodesy, 2008, 82(2): 65-82.

[23] AMIRI-SIMKOOEI A R. Application of least squares variance component estimation to errors-in-variables models[J]. Journal of Geodesy, 2013, 87(10): 935-944.

[24] XU Peiliang, LIU Jingnan, SHI Chuang. Total least squares adjustment in partial errors-in-variables models: algorithm and statistical analysis[J]. Journal of geodesy, 2012, 86(8): 661-675.

[25] YANG Y. Robust Estimation of Geodetic Datum Transformation[J]. Journal of geodesy, 1999, 73: 268-274.

[26] 王樂洋，許才軍.附有相對權(quán)比的總體最小二乘平差[J].武漢大學學報(信息科學報)，2011，36(8)：887-890．

[27] NERI F, SAITTA G, CHIOFALO S. An accurate and straightforward approach to line regression analysis of error-affected experimental data[J]. Journal of physics E scientific instruments, 1989, 8(22): 636-638.