陶葉青，蔡安寧，楊 娟

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

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

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

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

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

1 抗差總體最小二乘估計(jì)

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

(1)

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

(2)

(3)

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

(4)

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

2 基于方差分量估計(jì)的抗差算法

2.1 方差分量估計(jì)的最小二乘算法

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

Δ=Bx-L.

(5)

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

(6)

(7)

(8)

(9)

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

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

(10)

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

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

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

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

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

3 算例與分析

應(yīng)用直線擬合對(duì)建立的算法進(jìn)行驗(yàn)證，當(dāng)有n個(gè)觀測(cè)值時(shí)，將直線方程表示為

(11)

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

[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)南，等. 三維坐標(biāo)轉(zhuǎn)換的通用整體最小二乘算法[J].測(cè)繪學(xué)報(bào)，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] 王樂(lè)洋，許才軍，魯鐵定.邊長(zhǎng)變化反演應(yīng)變參數(shù)的總體最小二乘方法[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)報(bào))，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] 陳義，陸玨. 以三維坐標(biāo)轉(zhuǎn)換為例解算穩(wěn)健總體最小二乘方法[J].測(cè)繪學(xué)報(bào)，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] 龔循強(qiáng)，李志林. 穩(wěn)健加權(quán)總體最小二乘法[J]. 測(cè)繪學(xué)報(bào)，2014，43(9)：888-894.

[16] 龔循強(qiáng)，李志林.一種利用IGGⅡ方案的穩(wěn)健混合總體最小二乘方法[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)報(bào))，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]. 測(cè)繪學(xué)報(bào)，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ì)[J]. 測(cè)繪學(xué)報(bào)，2016，45(5)：552-559.

[21] 陶葉青，高井祥，姚一飛．基于中位數(shù)法的抗差總體最小二乘估計(jì)[J]．測(cè)繪學(xué)報(bào)，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] 王樂(lè)洋，許才軍.附有相對(duì)權(quán)比的總體最小二乘平差[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)報(bào))，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.