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

汪奇生 楊根新

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

?

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

汪奇生1楊根新2

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

基于懲罰函數和測量平差中權的思想，提出了附不等式約束的總體最小二乘平差模型，即利用懲罰函數對不等式約束方程構造約束權，通過零權和無限權將不等式約束轉換為等式約束，從而將不等式約束平差準則轉化為傳統的測量平差準則。同時，根據非線性最小二乘平差理論，用構造結構矩陣的方法來顧及系數矩陣的結構性，推導了附不等式約束的總體最小二乘迭代算法。該算法迭代格式與傳統的間接平差類似，只需經過若干次迭代便能得到最優解。

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

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

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

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

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

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

(1)

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

EIV模型的隨機模型為：

(2)

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

VA=DVa

(3)

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

(4)

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

(5)

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

(6)

這是一個附不等式約束的二次規劃問題。根據文獻[18]，按最優化理論，可以根據懲罰函數方法將上述問題轉化為無約束最優化問題：

(7)

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

(8)

式中，P0稱為約束權，其取值條件為：

(9)

根據式(8)、式(9)，當滿足約束條件時，P0(i)=0即為無效約束；當不滿足條件時，P0(i)=μ，為一個很大的數即為有效約束。通過平差中權的處理，使式(8)符合懲罰函數的特點。

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

(10)

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

(11)

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

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

展開后化簡得：

(12)

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

(13)

可將上式表示為：

(14)

根據間接平差原理，上式的解為：

(15)

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

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

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

3 實例分析

3.1 實例1

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

表1 算例原始數據

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

(16)

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

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

表2 估計結果

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

3.2 實例2

表3 線性回歸數據

表4 估計結果

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

4 結 語

本文根據最優化理論中懲罰函數方法，結合總體最小二乘平差準則提出的附有不等式約束的總體最小二乘迭代算法是可行的。該算法將不等式約束轉換為約束權的形式，迭代格式與間接平差相似，算法與傳統的平差方法接近，只需經過若干次迭代便可得到最優解。

本文的算法通過結構矩陣來處理系數矩陣含常數項或重復元素的情況，能較好地顧及到模型系數矩陣具有的結構特性。通過采用兩個算例進行分析，說明了本文算法的正確性，同時也驗證了本文算法能處理系數矩陣具有結構特性的問題。

[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]. 大地測量與地球動力學, 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].武漢大學學報:信息科學版, 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]馮光財,朱建軍.基于有效約束的附不等式約束平差的一種新算法[J].測繪學報, 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].測繪學報,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]. 測繪學報, 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]. 測繪學報,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

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

汪奇生，助理工程師，主要研究方向為大地測量數據處理,Email: wangqisheng0702@163.com。

10.14075/j.jgg.2016.12.015

1671-5942(2016)012-1100-05

P207

A