病態(tài)整體最小二乘的迭代正則化算法

孫同賀，羅志才，姚朝龍，宛家寬

?

病態(tài)整體最小二乘的迭代正則化算法

孫同賀1, 2，羅志才1, 3，姚朝龍1，宛家寬1

(1. 武漢大學(xué)測(cè)繪學(xué)院，武漢 430079；2. 內(nèi)蒙古科技大學(xué)礦業(yè)與煤炭學(xué)院，包頭014010；3. 武漢大學(xué)地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430079)

當(dāng)誤差含變量(EIV)模型的設(shè)計(jì)矩陣病態(tài)時(shí)，采用普通整體最小二乘(TLS)算法得不到穩(wěn)定的數(shù)值解。為了減弱病態(tài)性，在整體最小二乘準(zhǔn)則的基礎(chǔ)上附加解的二次范數(shù)約束，組成拉格朗日目標(biāo)函數(shù)，推導(dǎo)EIV模型的正則化整體最小二乘解(RTLS)。然后將RTLS的求解轉(zhuǎn)換為矩陣特征向量問(wèn)題，設(shè)計(jì)一個(gè)迭代方案逼近RTLS解。通過(guò)L曲線法求得正則化因子來(lái)確定正常數(shù)，從而避免人為選擇正常數(shù)的隨意性。數(shù)值實(shí)例表明，提出的迭代正則化算法是有效可行的。

EIV模型；病態(tài)問(wèn)題；正則化整體最小二乘；L曲線法；正常數(shù)

在信號(hào)處理、計(jì)算機(jī)視覺(jué)、圖像處理、通信工程以及大地測(cè)量與攝影測(cè)量等相關(guān)領(lǐng)域中，待求參數(shù)一般通過(guò)建立觀測(cè)模型來(lái)求解。當(dāng)觀測(cè)系統(tǒng)的設(shè)計(jì)矩陣和觀測(cè)量都含有誤差時(shí)，便形成了誤差含變量(Error-in-variables, EIV)模型。整體最小二乘(Total least squares, TLS)方法是解決EIV模型的基本途徑。TLS的算法主要包括奇異值分解類算法和拉格朗日迭代算法[1?2]，并擴(kuò)展到加權(quán)TLS[3]以及附有等式和不等式約束的情況[4?7]。在TLS解的統(tǒng)計(jì)性質(zhì)方面，主要研究了TLS解與LS解的關(guān)系[8]，系數(shù)陣的隨機(jī)誤差對(duì)加權(quán)LS解[9]的影響。對(duì)TLS在測(cè)繪領(lǐng)域的典型應(yīng)用也進(jìn)行了相關(guān)研究[10]。

當(dāng)設(shè)計(jì)陣的列向量間呈近似線性相關(guān)時(shí)，觀測(cè)量的微小變化會(huì)引起最小二乘解的劇烈振蕩，測(cè)量上一般采用截?cái)嗥娈愔祷蛘逿ikhonov正則化方法獲得穩(wěn)定的解[11?12]。根據(jù)奇異值分解，TLS解相較于LS解是一個(gè)降正則化過(guò)程，TLS解的條件數(shù)總是大于LS解的條件數(shù)[13]。GOLUB[14]提出了TLS的正則化解，給出了取不同的約束參數(shù)時(shí)，正則化整體最小二乘(Regularized total least squares, RTLS)解與TLS解的關(guān)系，并給出了一種直接算法。FIERRO等[15]給出了解RTLS的截?cái)嗥娈愔捣纸馑惴āCHAFFRIN等[16]基于Gauss-Helmert模型的混合近似解理論，提出了一種新的迭代計(jì)算方案，并運(yùn)用到考古學(xué)中來(lái)確定賽馬場(chǎng)中賽道的原點(diǎn)和半徑。BECK等[17]將RTLS問(wèn)題轉(zhuǎn)換成一個(gè)閉區(qū)間內(nèi)的單變量函數(shù)的最小化問(wèn)題。袁振超等[18]推導(dǎo)了等權(quán)條件下病態(tài)總體最小二乘的正則化解法。王樂(lè)詳?shù)萚19]采用嶺估計(jì)方法解決加權(quán)總體最小二乘平差的病態(tài)性問(wèn)題，給出了確定正則化參數(shù)的嶺跡法、廣義交叉核實(shí)法和L曲線法。又將參數(shù)作為虛擬觀測(cè)值，運(yùn)用廣義最小二乘準(zhǔn)則，推導(dǎo)了病態(tài)總體最小二乘問(wèn)題的虛擬觀測(cè)解法[20]。葛旭明等[21]提出了解RTLS的廣義正則化方法。然而，在應(yīng)用Tikhonov正則化方法解決病態(tài)問(wèn)題時(shí)，在該方法的約束條件中需人為引入一個(gè)正常數(shù)，但是正常數(shù)的引入并沒(méi)有絕對(duì)的標(biāo)準(zhǔn)，也沒(méi)有一個(gè)客觀的參考值，因而其值的引入存在較強(qiáng)的主觀隨意性。針對(duì)此問(wèn)題，引入了誤差約束的TLS正則化方法，在原有方程的基礎(chǔ)上加入2個(gè)不等式作為約束條件[22]。GUO等[23]直接給出了RTLS迭代算法的計(jì)算公式，而沒(méi)有給出詳細(xì)的推導(dǎo)過(guò)程。這些方法均是根據(jù)觀測(cè)量的精度及設(shè)計(jì)矩陣的形式人為選擇設(shè)計(jì)矩陣誤差上限值和正常數(shù)，因此具有隨機(jī)性。

本文作者基于文獻(xiàn)[23]中提出的RTLS的思想，詳細(xì)推導(dǎo)了病態(tài)EIV模型的正則化解及正常數(shù)與正則化因子的關(guān)系，然后將RTLS的求解轉(zhuǎn)換為求矩陣特征向量問(wèn)題，設(shè)計(jì)了一個(gè)迭代方案逼近RTLS解。通過(guò)L曲線法求得正則化因子來(lái)確定正常數(shù)，從而避免了人為選擇正常數(shù)的隨意性。最后用兩個(gè)算例說(shuō)明了本研究所提方法的可行性和有效性。

1 病態(tài)EIV模型的正則化方法

近年來(lái)，誤差含變量(EIV)模型在測(cè)量實(shí)際中廣泛出現(xiàn)。EIV模型采用奇異值分解算法求解時(shí)，發(fā)現(xiàn)它的整體最小二乘解相較于普通最小二乘解來(lái)說(shuō)，法方程主對(duì)角線元素分別減去增廣矩陣的最小奇異值，從而是一個(gè)降正則化(De-regularization)過(guò)程。王樂(lè)洋[13]證明了EIV模型的條件數(shù)要大于相應(yīng)的線性G-M模型。因此，有必要研究EIV模型病態(tài)時(shí)的數(shù)據(jù)處理問(wèn)題。

病態(tài)EIV模型可以表示為

在測(cè)量平差理論中，通過(guò)附加一些約束信息，可以削弱系統(tǒng)的病態(tài)性，保證估值的穩(wěn)定性。對(duì)于模型(1)所示的病態(tài)EIV模型，在觀測(cè)量和系數(shù)陣誤差殘差最小的基礎(chǔ)上，附加參數(shù)的二次范數(shù)約束，從而避免求解參數(shù)時(shí)因觀測(cè)量的微小變動(dòng)而造成解的巨大波動(dòng)。其平差準(zhǔn)則可表示成

式中：帽子符號(hào)表示估計(jì)值，波浪線符號(hào)表示預(yù)測(cè)值。根據(jù)式(4)和(5)可得

將式(9)、(10)代入式(8)，即

也可以表示成

整理上式得

將式(9)代入式(6)得

根據(jù)最優(yōu)化理論的Kuhn-Tucker條件和有效約束的概念，可以得到：若＞0，則，稱為有效約束；若，則＜為無(wú)效約束。由此，有

另外，由式(15)可以得到

令

則式(18)可以簡(jiǎn)寫成：

由式(6)易知：

將式(25)代入式(23)，有

將式(17)代入式(27)，可得

聯(lián)立式(21)和(29)，可得如下的矩陣形式：

式中：

式(21)的求解也可轉(zhuǎn)換為求矩陣特征值特征向量對(duì)的問(wèn)題，采用反冪法(Shifted inverse power)方法迭代求解，即求式(32)的特征向量：

定義第次迭代的殘差向量為

下面給出本研究所提方法的迭代步驟：

1) 給定約束陣，按照上述方法確定正則化因子，根據(jù)式(21)得到正則化整體最小二乘解；

2 算例分析

算例1：算例摘自文獻(xiàn)[24]，模型的設(shè)計(jì)矩陣和觀測(cè)值分別為

算例2：采用一組病態(tài)空間測(cè)邊網(wǎng)的算例，同時(shí)運(yùn)用不同方法進(jìn)行解算，進(jìn)而探討本方法在測(cè)量數(shù)據(jù)處理中的可行性。

表1 模型參數(shù)X的計(jì)算結(jié)果

模擬的空間測(cè)邊網(wǎng)算例，1、2、…9為9個(gè)已知點(diǎn)，其坐標(biāo)具體數(shù)據(jù)可參考文獻(xiàn)[22]。9個(gè)已知點(diǎn)到兩個(gè)未知點(diǎn)10、11的觀測(cè)距離也給出，假設(shè)兩個(gè)未知點(diǎn)真值坐標(biāo)為(0，0，0)和(7，10，?5)，它們之間的觀測(cè)距離為10,11=13.1078?m，各距離為等精度觀測(cè)，測(cè)距中誤差為±0.001?m。要求根據(jù)19個(gè)邊長(zhǎng)觀測(cè)值確定兩個(gè)未知點(diǎn)的坐標(biāo)。

經(jīng)計(jì)算，該測(cè)邊網(wǎng)觀測(cè)方程的系數(shù)矩陣嚴(yán)重病態(tài)，法方程的條件數(shù)為4.5886×103。在計(jì)算中，兩個(gè)未知點(diǎn)的坐標(biāo)近似值分別取為(0.01，?0.01，0.02)和(7.01，9.99，?5.01)。由文獻(xiàn)[22]可知，設(shè)計(jì)矩陣誤差上限可根據(jù)觀測(cè)量的精度確定，并在此基礎(chǔ)上略有加大。表2給出了各種算法的參數(shù)估計(jì)結(jié)果及與真值的差值范數(shù)，圖2給出了本方法得到的L曲線圖，橫縱坐標(biāo)軸同圖1。

通過(guò)公式推導(dǎo)及計(jì)算分析，可以得出如下結(jié)論：

2) 由表1和表2的統(tǒng)計(jì)結(jié)果可知，當(dāng)系數(shù)矩陣病態(tài)時(shí)，整體最小二乘解與真值的差值范數(shù)分別為6.7190和22.2787，較最小二乘解與真值的差值范數(shù)1.3088和10.4867要大，說(shuō)明整體最小二乘解受病態(tài)的影響更為嚴(yán)重；

3) 采用文獻(xiàn)[22]方法，設(shè)計(jì)矩陣誤差上限取不同值，得到的解與真值的差值范數(shù)分別為(0.8253，0.8257，0.8262)和(0.0925，0.0932，0.0938)；采用本研究方法得到的解與真值的差值范數(shù)分別為0.8234和0.0932。當(dāng)選擇較為合適的正常數(shù)時(shí)，文獻(xiàn)[22]中的方法和本研究方法解算結(jié)果相當(dāng)，但如果正常數(shù)選擇不合理，結(jié)果將有較大的偏差；

表2 模型參數(shù)X的計(jì)算結(jié)果

圖1 L曲線(算例1)

圖2 L曲線(算例2)

4) 在正則化最小二乘算法和正則化整體最小二乘算法中，通常假設(shè)正常數(shù)足夠小，計(jì)算過(guò)程中忽略了的確切值。本研究在迭代過(guò)程中考慮了的確切值。從解與真值的差值范數(shù)看，正則化整體最小二乘算法較之本研究方法得到的結(jié)果分別有0.0056和0.0038的偏差，這主要是因?yàn)檎?shù)帶來(lái)的影響；

5) 算例2中給定的初值并不理想，但利用本研究方法得到的解算結(jié)果要好于RLS和RTLS，由此說(shuō)明本研究方法對(duì)初值精度的要求不高，推廣了該方法的使用。

3 結(jié)論

1) 當(dāng)EIV模型設(shè)計(jì)矩陣病態(tài)時(shí)，普通的整體最小二乘算法得不到穩(wěn)定的解。在觀測(cè)量和系數(shù)陣誤差殘差最小的基礎(chǔ)上，附加參數(shù)的二次范數(shù)約束，推導(dǎo)了病態(tài)整體最小二乘的迭代正則化算法，以及正常數(shù)與正則化因子的等式關(guān)系，采用shifted inverse power算法，將RTLS的求解轉(zhuǎn)換為求矩陣的特征向量問(wèn)題。

2) 相比已有的正則化整體最小二乘算法，本研究提出的方法有兩個(gè)優(yōu)勢(shì)：首先，在迭代過(guò)程中考慮了的確切值；其次，通過(guò)L曲線法求得正則化因子來(lái)更合理地確定正常數(shù)，從而避免了不恰當(dāng)?shù)剡x取正常數(shù)，可能導(dǎo)致參數(shù)解存在較大偏差的情況。

REFERENCES：

[1] van HUFFEL S, VANDEWALLE J. The total least squares problem computational aspects and analysis[M]. Philadelphia: Society for Industrial and Applied Mathematics, 1991.

[2] GOLUB G H, van LOAN C F. An analysis of the total least squares problem[J]. SIAM Journal on numerical analysis, 1980, 17(6): 883?893.

[3] FANG Xing. Weighted total least squares solutions for applications in geodesy[D]. Hannover: Doctor’s Dissertation in Leibniz University, 2011.

[4] SCHAFFRIN B. A note on constrained total least-squares estimation[J]. Linear algebra and its applications, 2006, 417: 245?258.

[5] SCHAFFRIN B, FELUS Y A. An algorithmic approach to the total least-squares problem with linear and quadratic constraints[J]. Studia Geophysica et Geodaetica, 2009, 53: 1?16.

[6] ZHANG Song-lin, TONG Xiao-hua. A solution to EIV model with inequality constraints and its geodetic applications[J]. Journal of geodesy, 2013, 87: 23?28.

[7] FANG Xing. On non-combinatorial weighted total least squares with inequality constraints[J]. Journal of geodesy, 2014, 88: 805?816.

[8] van HUFFEL S, VANDEWALLE J. Algebraic connections between the least squares and total least squares problems[J]. Numerische Mathematik, 1989, 55: 431?449.

[9] XU Pei-liang, LIU Jing-nan. Effects of errors-in-variables on weighted least squares estimation[J]. Journal of Geodesy, 2014, 88: 705?716.

[10] 周擁軍, 鄧才華. 線性EIV模型的TLS估計(jì)及其典型應(yīng)用[J]. 中國(guó)有色金屬學(xué)報(bào), 2012, 22(3): 948?953. ZHOU Yong-jun, DENG Cai-hua. Extended total least squares problems for linear errors-in-variables model and its typical applications[J]. The Chinese Journal of Nonferrous Metals, 2012, 22(3): 948?953.

[11] XU Pei-liang. Truncated SVD methods for discrete linear ill-posed problems[J]. Geophysical Journal International, 1998, 135: 505?514.

[12] TIKHONOV A N. Solution of incorrectly formulated problems and the regularization method[J]. Soviet Math Dokl, 1963, 4: 1035?1038.

[13] 王樂(lè)洋. 總體最小二乘解性質(zhì)研究[J]. 大地測(cè)量與地球動(dòng)力學(xué), 2012, 32(5): 48?52. WANG Le-yang. Research on properties of total least squares estimation[J]. Journal of Geodesy and Geodynamics, 2012, 32(5): 48?52.

[14] GOLUB G H, HANSEN P C, O’LEARY D P. Tikhonov regularization and total least squares[J]. SIAM Journal on Matrix Analysis and Applications, 1999, 21(1): 185?194.

[15] FIERRO R D, GOLUB G H, HANSEN P C. Regularization by truncated total least squares[J]. SIAM Journal on Scientific Computing, 1997, 18(1): 1223?1241.

[16] SCHAFFRIN B, SNOW K. Total least-squares regularization of Tykhonov type and an ancient racetrack in Corinth[J]. Linear Algebra and Its Applications, 2010, 432: 2061?2076.

[17] BECK A, BEN-TAL A. On the solution of the Tikhonov regularization of the total least squares problem[J]. SIAM Journal on Optimization, 2006, 17(1): 98?118.

[18] 袁振超, 沈云中, 周澤波. 病態(tài)總體最小二乘模型的正則化算法[J]. 大地測(cè)量與地球動(dòng)力學(xué), 2009, 2(2): 131?134. YUAN Zhen-chao, SHEN Yun-zhong, ZHOU Ze-bo. Regularized total least squares solution to ill-posed error-in- variable model[J]. Journal of Geodesy and Geodynamics, 2009, 29(2): 131?134.

[19] 王樂(lè)洋, 許才軍, 魯鐵定. 病態(tài)加權(quán)總體最小二乘平差的嶺估計(jì)方法[J]. 武漢大學(xué)學(xué)報(bào)信息科學(xué)版, 2010, 35(11): 1346?1350. WANG Le-yang, XU Cai-jun, LU Tie-ding. Ridge estimation method in ill-posed weighted total least squares adjustment[J]. Geomatics and Information Science of Wuhan University, 2010, 35(11): 1346?1350.

[20] 王樂(lè)洋, 于冬冬. 病態(tài)總體最小二乘問(wèn)題的虛擬觀測(cè)解法[J]. 測(cè)繪學(xué)報(bào), 2014, 43(6): 575?581. WANG Le-yang, YU Dong-dong. Virtual observation method to ill-posed total least squares problem[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(6): 575?581.

[21] 葛旭明, 武吉倉(cāng). 病態(tài)總體最小二乘問(wèn)題的廣義正則化[J]. 測(cè)繪學(xué)報(bào), 2012, 41(3): 372?377. GE Xu-ming, WU Ji-cang. Generalized regularization to ill-posed total least squares problem[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(3): 372?377.

[22] 葛旭明, 武吉倉(cāng). 誤差限的病態(tài)總體最小二乘解算[J]. 測(cè)繪學(xué)報(bào), 2013, 42(2): 196?202. GE Xu-ming, WU Ji-cang. A regularization method to ill-posed total least squares with error limits[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(2): 192?202.

[23] GUO Hong-bin, RENAUT R A. A regularized total least squares algorithm[C]// Total Least Squares and Errors-in-Variables Modeling. Copenhagen: Kluwer Academic Publishers, 2002: 57?66.

[24] 魯鐵定, 寧津生.總體最小二乘平差理論及其應(yīng)用[M]. 北京: 中國(guó)科學(xué)技術(shù)出版社, 2011. LU Tie-ding, Ning Jin-sheng. Total least squares adjustment theory and its applications[M]. Beijing: China Science and Technology Press, 2011.

(編輯 王 超)

Iterative and regularized algorithm to ill-posed total least squares

SUN Tong-he1, 2, LUO Zhi-cai1, 3, YAO Chao-long1, YUAN Jia-kuan1

(1. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China; 2. Mining and Coal Institute, Inner Mongolia University of Science and Technology, Baotou 014010, China; 3. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China)

When the design matrix of errors-in-variables (EIV) model was ill-conditioned, the ordinary total least squares (TLS) solution was unstable. In order to weaken the ill-conditioning, an Euclid norm constraint of the solution was added to the TLS minimization rule. Then, the Lagrange objective function was formed and the regularized total least squares (RTLS) solution was deduced. Afterwards, the RTLS was transformed to a problem of looking for a matrix’s eigenvector. An iterative program was designed to approximate the solution. The L-curve method was used to choose the regularization factor to determine the positive constant, which can avoid the subjective decision. The simulations show the efficiency and feasibility of the algorithm.

EIV model; ill-posed problem; regularized total least squares; L-curve method; positive constant

Project(41474006) supported by the National Natural Science Foundation of China

2015-12-04; Accepted date:2016-04-10

SUN Tong-he; Tel: +86-472-5951692; E-mail: suntonghe202@163.com

1004-0609(2016)-10-2174-07

P207

A

國(guó)家自然科學(xué)基金資助項(xiàng)目(41474006)

2015-12-04；

2016-04-10

孫同賀，講師，博士；電話：0472-5951692；E-mail: suntonghe202@163.com