附不等式約束的總體最小二乘迭代算法

汪奇生 楊根新

1 湖南軟件職業(yè)學(xué)院建筑工程學(xué)院，湘潭市開源路1號，411100 2 云南國土資源職業(yè)學(xué)院測繪地理信息學(xué)院，昆明市經(jīng)牛路2號，650217

?

附不等式約束的總體最小二乘迭代算法

汪奇生1楊根新2

1 湖南軟件職業(yè)學(xué)院建筑工程學(xué)院，湘潭市開源路1號，411100 2 云南國土資源職業(yè)學(xué)院測繪地理信息學(xué)院，昆明市經(jīng)牛路2號，650217

基于懲罰函數(shù)和測量平差中權(quán)的思想，提出了附不等式約束的總體最小二乘平差模型，即利用懲罰函數(shù)對不等式約束方程構(gòu)造約束權(quán)，通過零權(quán)和無限權(quán)將不等式約束轉(zhuǎn)換為等式約束，從而將不等式約束平差準(zhǔn)則轉(zhuǎn)化為傳統(tǒng)的測量平差準(zhǔn)則。同時(shí)，根據(jù)非線性最小二乘平差理論，用構(gòu)造結(jié)構(gòu)矩陣的方法來顧及系數(shù)矩陣的結(jié)構(gòu)性，推導(dǎo)了附不等式約束的總體最小二乘迭代算法。該算法迭代格式與傳統(tǒng)的間接平差類似，只需經(jīng)過若干次迭代便能得到最優(yōu)解。

不等式約束；EIV模型；總體最小二乘；迭代算法；懲罰函數(shù)

總體最小二乘(total least squares)是一種能同時(shí)考慮系數(shù)矩陣誤差的方法[1]，受到各領(lǐng)域?qū)W者的廣泛關(guān)注。在測量數(shù)據(jù)處理中，總體最小二乘估計(jì)方法對應(yīng)的平差模型為EIV(errors in variables)模型。對于EIV模型的解算，國內(nèi)外學(xué)者進(jìn)行了深入研究[2-9]。其中，文獻(xiàn)[2]運(yùn)用拉格朗日原理首次提出總體最小二乘的迭代法，文獻(xiàn)[3]針對線性回歸系數(shù)矩陣含有常數(shù)列提出了其總體最小二乘解，文獻(xiàn)[4-8]研究了加權(quán)總體最小二乘算法并應(yīng)用于測量數(shù)據(jù)處理。除此之外，一些學(xué)者還研究了擴(kuò)展總體最小二乘的一些其他算法[9]。

以上算法都沒有考慮參數(shù)估計(jì)時(shí)的先驗(yàn)信息。當(dāng)存在某些先驗(yàn)信息時(shí)，可根據(jù)先驗(yàn)信息對參數(shù)附加某種約束。如果約束是等式，則可以構(gòu)建附有等式約束的總體最小二乘模型(equality constrained EIV,ECEIV)。對于ECEIV模型的解算，文獻(xiàn)[10-12]進(jìn)行了詳細(xì)論述。如果約束是不等式，則可以構(gòu)建附有不等式約束的總體最小二乘模型(inequality constrained EIV,ICEIV)。對于附有不等式約束的平差問題，基于最小二乘的研究成果較多[13-18]，而基于總體最小二乘的研究則較少[19-22]。其中，文獻(xiàn)[19]采用線性互補(bǔ)方法通過排列組合進(jìn)行求解，首次提出了附不等式約束的總體最小二乘解算方法；文獻(xiàn)[20]采用窮舉法，將不等式約束轉(zhuǎn)換為等式約束問題進(jìn)行求解。這兩種方法算法效率低，且僅適用于系數(shù)矩陣元素全部為隨機(jī)元素的EIV 模型以及等精度、不相關(guān)觀測值[21]。因此，文獻(xiàn)[21]采用懲罰函數(shù)的方法提出一種適應(yīng)于系數(shù)矩陣為任意形式的算法；文獻(xiàn)[22]則采用線性互補(bǔ)方法提出一種計(jì)算效率高的迭代算法。盡管上述文獻(xiàn)分別提出了一種附有不等式約束的總體最小二乘算法，而且算法效率從低到高，算法從僅適應(yīng)于系數(shù)矩陣為一般形式到任意形式，但都存在算法復(fù)雜、計(jì)算困難的問題，而且這些算法都不是基于傳統(tǒng)的平差方法，不利于測量人員理解。

本文根據(jù)文獻(xiàn)[18]處理最小二乘的附不等式約束平差時(shí)的思想，基于懲罰函數(shù)處理不等式約束的方法，結(jié)合傳統(tǒng)測量平差中權(quán)的思想，提出一種類似間接平差模型的簡單迭代算法。與文獻(xiàn)[21]采用的懲罰函數(shù)方法不同的是，本文將懲罰函數(shù)處理不等式約束問題的思想與測量平差方法相結(jié)合。該算法利用構(gòu)造結(jié)構(gòu)矩陣來顧及系數(shù)矩陣的結(jié)構(gòu)性，能處理系數(shù)矩陣為任意形式的ICEIV模型，算法迭代格式簡單易懂，易于測量人員理解和掌握。

1 附不等式約束的總體最小二乘平差模型(ICEIV)

附不等式約束的總體最小二乘平差模型(ICEIV模型)為[20]：

(1)

式中，L為m×1的觀測向量，VL為m×1的觀測值改正向量，A為m×n的系數(shù)矩陣，EA為m×n的系數(shù)改正矩陣，X為n×1的參數(shù)估值向量，G為k×n不等式約束方程的系數(shù)矩陣且為行滿秩矩陣，W為k×1的常數(shù)向量。

EIV模型的隨機(jī)模型為：

(2)

然而在實(shí)際處理時(shí)，EIV模型的系數(shù)矩陣并不是全部都含有誤差，因此本文提出用結(jié)構(gòu)矩陣的方法來顧及這一問題，用結(jié)構(gòu)矩陣來表示系數(shù)矩陣的改正向量VA。設(shè)VA中有t×1的元素含有誤差，則：

VA=DVa

(3)

式中，Va是t×1的系數(shù)矩陣中含誤差的非重復(fù)元素的改正數(shù)(數(shù)量為t個(gè))，D是mn×t的結(jié)構(gòu)矩陣。因此，ICEIV平差模型可以表示為：

(4)

式中，vec-1(·)表示vec(·)的逆運(yùn)算。ICEIV的隨機(jī)模型為：

(5)

按照總體最小二乘原理，ICEIV模型的平差準(zhǔn)則為：

(6)

這是一個(gè)附不等式約束的二次規(guī)劃問題。根據(jù)文獻(xiàn)[18]，按最優(yōu)化理論，可以根據(jù)懲罰函數(shù)方法將上述問題轉(zhuǎn)化為無約束最優(yōu)化問題：

(7)

式中，P(x)為懲罰函數(shù)。式(7)的基本思想是對目標(biāo)函數(shù)式(6)中不滿足約束條件的點(diǎn)給予懲罰，即給P(x)取一個(gè)很大的值；當(dāng)滿足約束條件時(shí)，則P(x)取值為0，轉(zhuǎn)化為無約束的普通總體最小二乘問題。在最優(yōu)化計(jì)算中，懲罰函數(shù)具體形式一般憑經(jīng)驗(yàn)確定，也可由先驗(yàn)信息的置信度決定。對于式(6)的不等式約束，令V0=GX-W，則可構(gòu)造懲罰函數(shù)：

(8)

式中，P0稱為約束權(quán)，其取值條件為：

(9)

根據(jù)式(8)、式(9)，當(dāng)滿足約束條件時(shí)，P0(i)=0即為無效約束；當(dāng)不滿足條件時(shí)，P0(i)=μ，為一個(gè)很大的數(shù)即為有效約束。通過平差中權(quán)的處理，使式(8)符合懲罰函數(shù)的特點(diǎn)。

由此，附不等式約束的總體最小二乘平差準(zhǔn)則可由式(6)轉(zhuǎn)化為：

(10)

通過上文分析，實(shí)際上附不等式約束的總體最小二乘平差模型可以轉(zhuǎn)化為等式約束的總體最小二乘平差模型，即式(4)可以等價(jià)為：

(11)

在式(10)的平差準(zhǔn)則下，通過迭代方法可以求得式(11)的最優(yōu)解。

2 附不等式約束的總體最小二乘迭代算法(ICTLS)

展開后化簡得：

(12)

式中，?為矩陣的可內(nèi)克積。根據(jù)總體最小二乘原理，可將式(12)寫成間接平差形式：

(13)

可將上式表示為：

(14)

根據(jù)間接平差原理，上式的解為：

(15)

式(14)即為普通的最小二乘間接平差模型。要求附不等式約束的總體最小二乘解，需要通過式(15)進(jìn)行迭代計(jì)算。式(13)、(14)即為迭代的基本格式。

本文的附不等式約束的總體最小二乘迭代算法的具體解算步驟為：

4)重復(fù)第3步，直到兩次計(jì)算的參數(shù)之差小于給定的限差，則停止迭代，輸出參數(shù)值。

3 實(shí)例分析

3.1 實(shí)例1

實(shí)例數(shù)據(jù)選自文獻(xiàn)[13]，如表1所示。表中的數(shù)據(jù)包括系數(shù)矩陣和常數(shù)向量的值，并附有3個(gè)不等式約束，同時(shí)對xi的取值有限制。

表1 算例原始數(shù)據(jù)

在表1中，共有11個(gè)不等式約束，其中對xi的取值限制可以表示為8個(gè)不等式約束，對應(yīng)的x1約束可以表示成：

(16)

根據(jù)式(16)，可以類似將x2、x2、x3的取值限制表示成不等式約束。

分別采用總體最小二乘法、文獻(xiàn)[20]法、文獻(xiàn)[21]法、文獻(xiàn)[22]法以及本文的方法進(jìn)行參數(shù)求解，其中采用本文方法構(gòu)造的結(jié)構(gòu)矩陣為單位矩陣D=eye(20)，取迭代停止條件ε=10-10。參數(shù)求解的結(jié)果如表2所示。

表2 估計(jì)結(jié)果

從表2可以發(fā)現(xiàn)，由于本例的系數(shù)矩陣全部含有誤差，因此采用本文方法與其他3種方法得到的參數(shù)估值是完全相同的，即滿足先驗(yàn)信息的約束條件，而總體最小二乘解則不能滿足約束條件,這說明本文所述的附有不等式約束的總體最小二乘迭代算法是可行的。相比之下，本文采用的迭代算法既能彌補(bǔ)文獻(xiàn)算法計(jì)算量大的缺陷，又能通過結(jié)構(gòu)矩陣來顧及系數(shù)矩陣的結(jié)構(gòu)性，同時(shí)算法的迭代形式簡單，易于編程實(shí)現(xiàn)。

3.2 實(shí)例2

表3 線性回歸數(shù)據(jù)

表4 估計(jì)結(jié)果

從表4可以看出，對于系數(shù)矩陣具有結(jié)構(gòu)性的附不等式約束的總體最小二乘模型，本文算法求得的參數(shù)估值與文獻(xiàn)[21]、文獻(xiàn)[22]的結(jié)果相同，都滿足約束條件，而總體最小二乘解則不滿足。這再次驗(yàn)證了本文算法可以解算系數(shù)矩陣具有任意結(jié)構(gòu)特性的附不等式約束的總體最小二乘模型。

4 結(jié) 語

本文根據(jù)最優(yōu)化理論中懲罰函數(shù)方法，結(jié)合總體最小二乘平差準(zhǔn)則提出的附有不等式約束的總體最小二乘迭代算法是可行的。該算法將不等式約束轉(zhuǎn)換為約束權(quán)的形式，迭代格式與間接平差相似，算法與傳統(tǒng)的平差方法接近，只需經(jīng)過若干次迭代便可得到最優(yōu)解。

本文的算法通過結(jié)構(gòu)矩陣來處理系數(shù)矩陣含常數(shù)項(xiàng)或重復(fù)元素的情況，能較好地顧及到模型系數(shù)矩陣具有的結(jié)構(gòu)特性。通過采用兩個(gè)算例進(jìn)行分析，說明了本文算法的正確性，同時(shí)也驗(yàn)證了本文算法能處理系數(shù)矩陣具有結(jié)構(gòu)特性的問題。

[1] Golub G H,Vanloan C F．An Analysis of the Total Least Squares Problem[J]．SIAM Journal on Numerical Analysis,1980,17(6):883G893

[2] Schaffrin B. A Note on Constrained Total Least-Squares Estimation[J]. Linear Algebra and Its Applications, 2006, 417(1): 245-258

[3] 汪奇生，楊德宏，楊建文. 基于總體最小二乘的線性回歸迭代算法[J]. 大地測量與地球動(dòng)力學(xué), 2013, 33(6): 112-114 (Wang Qisheng，Yang Dehong，Yang Jianwen. An Iteration Algorithm of Linear Regression Based on Total Least Squares[J]. Journal of Geodesy and Geodynamics，2013, 33(6): 112-114)

[4] Schaffrin B,Wieser A. On Weighted Total Least-Squares Adjustment for Linear Regression[J].Journal of Geodesy,2008,82(7): 415-421

[5] Shen Y Z,Li B F,Chen Y. An Iterative Solution of Weighted Total Least-Squares Adjustment[J]. Journal of Geodesy,2010, 85(4): 229-238

[6] Mahboub V.On Weighted Total Least-Squares for Geodetic Transformations[J].Journal of Geodesy, 2012,86(5): 359-367[7] 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

[8] Xu P L,Liu J N,Shi C.Total Least Squares Adjustment in Partial Errors-in-Variables Models:Algorithm and Statistical Analysis[J]. Journal of Geodesy,2012, 86(8):661-675

[9] 劉經(jīng)南,曾文憲,徐培亮.整體最小二乘估計(jì)的研究進(jìn)展[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版, 2013, 38(5): 505-512(Liu Jingnan, Zeng Wenxian, Xu Peiliang. Overviwe of Total Least Squares Methods[J]. Geomatics and Information Science of WuhanUniversity,2013,38(5):505-512)

[10]Schaffrin B,Felus Y A．On Total Least Squares Adjustment with Constraints[C].International Associationof Geodesy Symposia,Sapporo,2005

[11]Mahboub V, Sharifi M A. On Weighted Total Least-Squares with Linear and Quadratic Constraints[J]. Journal of Geodesy, 2013, 87(3): 279-286

[12]Fang X. A Structured and Constrained Total Least-Squares Solution with Cross-Covariances[J]. Studia Geophysica et Geodaetica, 2013: 1-16

[13]Peng J，Guo C, Zhan H． An Aggregate Constraint Method for Inequality-Constrained Least Squares Problem[J].Journal of Geodesy，2006，79(12): 705-713

[14]Zhu J,Santerre R,Chang X W.A Bayesian Method for Linear Inequality Constrained Adjustment and Its Application to GPS Positioning[J].Journal of Geodesy,2005,78(9):528-534

[15]馮光財(cái),朱建軍.基于有效約束的附不等式約束平差的一種新算法[J].測繪學(xué)報(bào), 2007,36(2):119-122(Feng Guangcai,Zhu Jianjun.A New Approach to Inequality Constrained Least-Squares Adjustment[J].Acta Geodaeticaet Cartographica Sinica,2007,36(2):119-122)

[16]朱建軍,歐陽文森,文小岳.基于遺傳算法解決附有不等式約束的最小二乘平差問題的研究[J].工程勘察,2006(3):61-64(Zhu Jianjun,Ouyang Wensen,Wen Xiaoyue.Solving the LICA Problem by GA Method [J].Journal of Geotechnical Investigation & Surveying,2006(3):61-64)

[17]宋迎春,左廷英,朱建軍.帶有線性不等式約束平差模型的算法研究[J].測繪學(xué)報(bào),2008,37(4):433-437(Song Yingchun,Zuo Tingying,Zhu Jianjun.Research on Algorithm of Adjustment Model with Linear Inequality Constrained Parameters[J]. Acta Geodaeticaet Cartographica Sinica,2008,37(4):433-437)

[18]朱建軍,謝建. 附不等式約束平差的一種簡單迭代算法[J]. 測繪學(xué)報(bào), 2011, 40(2): 209-211(Zhu Jianjun,Xie Jian. A Simple Iterative Algorithm for Inequality Constrained Adjustment[J]. Acta Geodaeticaet Cartographica Sinica, 2011, 40(2): 209-211)

[19]Moor B.Total Linear Least Squares with Inequality Constraints[R].Department of Electrical Engineering,Delft University,1990

[20]Zhang S L, Tong X H, Zhang K．A Solution to EIV Model with Inequality Constraints and Its Geodetic Applications[J]. Journal of Geodesy,2013, 87(1):23-28

[21]曾文憲,方興,劉經(jīng)南,姚宜斌.附有不等式約束的加權(quán)整體最小二乘算法[J]. 測繪學(xué)報(bào),2014,43(10):1 013-1 018 (Zeng Wenxian,Fang Xing,Liu Jingnan,et al．Weighted Total Least Squares Algorithm with Inequality Constraints[J]．Acta Geodaetica et Cartographica Sinica,2014,43(10): 1 013-1 018)

[22]Zeng W X,Liu J N,Yao Y B. On Partial Errors-in-Variables Models with Inequality Constraints of Parameters and Variables[J].Journal of Geodesy，2015，89(2): 111-119

About the first author:WANG Qisheng，assistant engineer，majors in geodetic date processing,E-mail:wangqisheng0702@163.com.

1 Iteration Algorithm of Total Least Squares with Inequality Constraints

WANGQisheng1YANGGenxin2

1 Institute of Civil Engineering, Hunan Software Vocational Institute, 1 Kaiyuan Road, Xiangtan 411100, China 2 College of Geomatics and Geoinformation, Yunnan Land and Resources Vocational College, 2 Jingniu Road,Kunming 650217,China

Based on penalty function and weight of adjustment, an inequality constraints EIV (ICEIV) model is presented. The model utilizes penalty function to construct the constraint weight for the constraint equations and transforms the inequality constraint into the equality constraint by the zero or infinite weight. So, it can transform the inequality constraint adjustment criteria into the classical adjustment criteria. Therefore, a new iteration algorithm of total least squares with inequality constraints is deduced by the nonlinear least squares adjustment theory; the method uses a structured matrix to consider the repetitive elements and constant terms.

inequality constraints；errors-in-variables model；total least squares； iteration algorithm； penalty function

Scientific Research Foundation of Education Bureau of Hunan Province,No.15C0741;Scientific Research Foundation of Education Bureau of Yunnan Province,No.2016ZZX252.

2015-12-26

項(xiàng)目來源：湖南省教育廳科研項(xiàng)目(15C0741);云南省教育廳科研基金(2016ZZX252)。

汪奇生，助理工程師，主要研究方向?yàn)榇蟮販y量數(shù)據(jù)處理,Email: wangqisheng0702@163.com。

10.14075/j.jgg.2016.12.015

1671-5942(2016)012-1100-05

P207

A