附有不等式約束的加權(quán)整體最小二乘算法

曾文憲,方 興,劉經(jīng)南,2,姚宜斌

1.武漢大學(xué)測繪學(xué)院,湖北武漢 430079;2.武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心,湖北武漢 430079

附有不等式約束的加權(quán)整體最小二乘算法

曾文憲1,方 興1,劉經(jīng)南1,2,姚宜斌1

1.武漢大學(xué)測繪學(xué)院,湖北武漢 430079;2.武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心,湖北武漢 430079

針對現(xiàn)有附有不等式約束的整體最小二乘算法的缺陷,以partial EIV(errors-in-variables)模型為基礎(chǔ),在整體最小二乘準(zhǔn)則下,通過將附有不等式約束的EIV模型的求解轉(zhuǎn)換為標(biāo)準(zhǔn)的附有不等式約束的最優(yōu)化問題,并采用懲罰函數(shù)法等方法得到了附有不等式約束的加權(quán)整體最小二乘新算法。新算法將現(xiàn)有算法的特殊權(quán)陣限制條件擴展到了一般性權(quán)矩陣,將要求系數(shù)矩陣元素全部隨機的限定條件擴展到了可同時包含隨機和非隨機元素的一般情況,并且新算法解決了現(xiàn)有算法計算量受制于約束方程數(shù)量的缺陷。實例計算表明,本文提出的算法簡單、有效,具有普遍適用性。

整體最小二乘估計;EIV模型;不等式約束;非線性算法

1 引 言

整體最小二乘估計(total least squares, TLS)作為EIV[1](errors-in-variables)模型的嚴(yán)密估計方法,目前已廣泛應(yīng)用于大地測量等眾多科學(xué)研究和工程應(yīng)用領(lǐng)域。文獻(xiàn)[2]提出了整體最小二乘準(zhǔn)則。文獻(xiàn)[3]基于正交回歸原理推導(dǎo)了TLS數(shù)值算法。假定觀測值不相關(guān)情況下,文獻(xiàn)[4]首次提出了真正統(tǒng)計意義上的TLS算法。文獻(xiàn)[5]提出了著名的奇異值分解(singular value decomposition,SVD)算法。文獻(xiàn)[6]研究了穩(wěn)健整體最小二乘算法(robust TLS)等。文獻(xiàn)[7]證明了當(dāng)觀測數(shù)趨于無窮大時,TLS估計具有漸進(jìn)無偏性。大地測量領(lǐng)域針對普遍存在的觀測值不等精度、相關(guān)的情況,文獻(xiàn)[8—11]提出了加權(quán)TLS算法(weighted TLS,WTLS),其中,文獻(xiàn)[11]研究了最一般性權(quán)矩陣條件下的WTLS算法。文獻(xiàn)[1]提出了基于partial EIV模型(PEIV)的WTLS算法,能夠?qū)⒔Y(jié)構(gòu)性等各種形式的系數(shù)矩陣納入統(tǒng)一的模型形式求解。其他TLS算法的發(fā)展情況見文獻(xiàn)[12]。

當(dāng)參數(shù)估計存在先驗信息時,EIV模型擴展為附有約束的EIV模型。相當(dāng)多的文獻(xiàn)研究了附有等式約束的EIV模型(equality-constrained EIV,ECEIV),如文獻(xiàn)[13]、文獻(xiàn)[14]以及文獻(xiàn)[15]分別提出了3種不同的附有等式約束的TLS算法(equality-constrained TLS,ECTLS)。

某些先驗信息要求對模型附加不等式約束條件,就形成了附有不等式約束的EIV模型(inequality-constrained EIV,ICEIV),如變形監(jiān)測中某因素引起的變形系數(shù)理論上要大于某一數(shù)值,應(yīng)對該估計值進(jìn)行最小值約束;坐標(biāo)轉(zhuǎn)換模型的TLS算法中要求強制附合到某中心點,即該點坐標(biāo)的改正數(shù)不能超出一定數(shù)值范圍等。關(guān)于附有不等式約束的EIV模型目前僅檢索到兩篇研究文獻(xiàn),文獻(xiàn)[16]運用拉格朗日乘數(shù)法將ICTLS問題轉(zhuǎn)化為廣義線性互補問題(linear complementarity problem,LCP),再通過排列組合的方法進(jìn)行求解。文獻(xiàn)[17]提出了基于窮舉法的附有不等式約束的整體最小二乘算法(inequality-constrained TLS,ICTLS)。上述文獻(xiàn)開啟了ICTLS算法研究的先河,但存在以下問題有待解決:①算法的計算量隨約束方程個數(shù)的增長呈指數(shù)增長,如當(dāng)約束方程個數(shù)為50時,組合或搜索次數(shù)達(dá)到約250≈1015,因此,約束方程較多時,算法的計算量急劇增長甚至導(dǎo)致無法計算;②算法僅適用于系數(shù)矩陣元素全部為隨機元素的EIV模型以及等精度、不相關(guān)觀測值。

本文以partial EIV模型為基礎(chǔ),研究了附有不等式約束的加權(quán)整體最小二乘(inequality-constrained weighted TLS,ICWTLS)算法。本文提出的ICWTLS算法適用于隨機和非隨機元素并存、或呈現(xiàn)結(jié)構(gòu)性特征的一般性系數(shù)矩陣以及不等精度、相關(guān)觀測值,并且算法的計算效率遠(yuǎn)高于基于窮舉法的ICTLS算法,計算量不受約束方程個數(shù)的制約。

2 附有不等式約束的partial-EIV模型

ICEIV模型的數(shù)學(xué)表達(dá)式為[7,17])

式中,y表示n×1的觀測向量;ey表示y的n×1隨機誤差向量,且ey期望為0、協(xié)因數(shù)陣Qy＝w－1σ2(w和σ2分別表示觀測向量的權(quán)陣以及單位權(quán)方差);β表示t×1的參數(shù)向量;A表示n×t的系數(shù)矩陣;EA表示A的n×t隨機誤差矩陣, EA的期望為0且協(xié)因數(shù)陣為QA＝ω－1σ2(ω表示系數(shù)矩陣觀測元素的權(quán)陣),通常假定EA與ey不相關(guān);G表示不等式約束方程的k×t系數(shù)矩陣;z表示k×1的常數(shù)向量。

現(xiàn)有ICTLS算法[16-17]假定式(1)中觀測數(shù)據(jù)的權(quán)陣ω和w均為單位陣,并且A中元素全部為隨機量。當(dāng)A中存在非隨機的固定元素或者呈現(xiàn)結(jié)構(gòu)性特征時,必須對系數(shù)矩陣的協(xié)因數(shù)陣[18]或者對系數(shù)矩陣[15]進(jìn)行特殊處理后求解。為了得到一般性系數(shù)矩陣下的統(tǒng)一算法,以partial EIV模型[1]為基礎(chǔ),將傳統(tǒng)的ICEIV模型(1)改寫為如下形式

式(2)構(gòu)成了附有不等式約束的partial EIV模型(inequality-constvained partial EIV,ICPEIV),與傳統(tǒng)的ICEIV模型[16-17]比較,其優(yōu)勢主要體現(xiàn)為:①將ICEIV模型中系數(shù)矩陣元素全部為隨機量的限定擴展到了同時包含隨機和非隨機元素的一般情況;②ICPEIV模型可解析結(jié)構(gòu)性系數(shù)矩陣,保證了算法的統(tǒng)一性;③ICPEIV模型A中參與平差計算的待估量個數(shù)m小于等于ICEIV中A中元素個數(shù)nt,尤其當(dāng)系數(shù)矩陣中的隨機量個數(shù)較少時,ICPEIV模型形式可以大大減少待估計量。因此,ICPEIV模型更具一般性,以下在不限定觀測數(shù)據(jù)權(quán)陣的等精度和相關(guān)性的一般情況下,即ICPEIV模型(2)中權(quán)矩陣ω和w為任意正定對稱陣,筆者提出了普遍適用的附有不等式約束的加權(quán)整體最小二乘新算法。

3 (ICWTLS)算法

可以求出其最優(yōu)估計值。由式(2)可得

將ey和ea代入式(3),則ICPEIV模型(2)的ICWTLS算法轉(zhuǎn)化為以下附有不等式約束的最優(yōu)化問題

通過將ICPEIV模型的估計轉(zhuǎn)換為附有不等式約束的最優(yōu)化問題,避免了現(xiàn)有算法采用窮舉法引起的算法受限于不等式方程個數(shù)的缺陷。根據(jù)最優(yōu)化理論,提出了基于罰函數(shù)法[19]的ICWTLS新算法。

式中,μ表示懲罰參數(shù)(μ＞0);ci(β)＝max[0,－Giβ＋zi](i＝1,2,…,k);Gi為G的第i行;zi為z的第i個元素。式(5)存在一階導(dǎo)數(shù)且一階導(dǎo)數(shù)連續(xù)。

基于罰函數(shù)法的ICWTLS算法計算過程如下。

(3)如果參數(shù)解滿足收斂條件,結(jié)束計算,轉(zhuǎn)至步驟(4)。

(4)按給定規(guī)則增大μ0,轉(zhuǎn)至步驟(2)進(jìn)行迭代計算。

對于ICPEIV模型估計轉(zhuǎn)換后的最優(yōu)化式(4),除上述基于罰函數(shù)的ICWTLS算法外,同樣可以采用最優(yōu)估計理論的有效集法、序列二次規(guī)劃法、內(nèi)點算法等[19-20]得到相應(yīng)的ICWTLS算法?；谏鲜霾煌顑?yōu)估計方法的ICWTLS算法的特點是下一步要討論的問題。

4 實例分析

為了說明本文提出的附有不等式約束的整體最小二乘算法的應(yīng)用,筆者共模擬兩個實例進(jìn)行了計算。實例1的主要目的是驗證和比較ICWTLS新算法與現(xiàn)有基于窮舉法的ICTLS算法[17]的優(yōu)缺點及計算效率,因此,筆者設(shè)計了一組觀測數(shù)據(jù)等精度并且系數(shù)矩陣全部為隨機元素的ICEIV模型數(shù)據(jù)。實例2選擇了附有不等式約束的平面擬合模型,該模型的系數(shù)矩陣同時存在隨機和非隨機元素,現(xiàn)有算法無法解算,通過模擬一組不等精度的試驗數(shù)據(jù),采用ICWTLS算法求出了模型的參數(shù)解。

4.1 實例1(等精度、系數(shù)矩陣為隨機元素)

實例1的模擬數(shù)據(jù)見表1,模型包含7×1待估參數(shù)向量β以及70×1待估系數(shù)向量ˉa,系數(shù)矩陣A的全部元素為隨機量,A和y中所有觀測元素不相關(guān)且中誤差均為0.1,其中參數(shù)真值β～見表2。為了改進(jìn)模型的估計結(jié)果,利用參數(shù)的先驗信息設(shè)計了18個不等式約束方程(見表1),其中,0≤βi≤0.8(i＝1,2,3,4)以及－0.5≤βj≤0 (j＝5,6,7)可分別表示為

式中,I4和I3分別表示4階和3階的單位陣。

表1 附有不等式約束的平差模型數(shù)據(jù)Tab.1 Data set of the inequality-constrained adjustment model

(2)從算法的計算量比較,本文算法只需迭代計算23次,而窮舉法理論上排列組合所需計算次數(shù)為218≈260 000。當(dāng)模型約束方程個數(shù)較多的情況下,基于窮舉法的ICTLS算法計算量要遠(yuǎn)大于本文提出的ICWTLS算法。當(dāng)不等式約束方程過多時,基于窮舉法的算法甚至無法在有效時間內(nèi)求解。

表2 ICWTLS算法估計結(jié)果Tab.2 Estimation of ICWTLS

4.2 實例2(附有不等式約束的平面擬合模型)

線性回歸模型是測繪領(lǐng)域常用的基本EIV模型之一,實例2選用平面擬合模型說明ICWTLS算法的應(yīng)用,平面擬合模型式(1)中觀測向量、系數(shù)矩陣、參數(shù)向量等形式如下

由上式可以看到平面擬合模型的系數(shù)矩陣同時包含了非隨機和隨機元素,即第1列元素為已知量,第2列和第3列是由觀測數(shù)據(jù)構(gòu)成的隨機量,現(xiàn)有附有不等式約束的整體最小二乘算法無法處理這類情況。本文算法可以解算任意結(jié)構(gòu)性系數(shù)矩陣形式的EIV模型,將上式表示為partial EIV模型式(2)的形式,式中的矩陣和向量可表示為

式中,1n×1表示元素均為1的n×1單位向量; 02n×1表示元素均為0的2n×1向量;0n×n表示n維的零方陣;I2n×2n表示2n維的單位陣;bn×1＝ [b1b2…bn]T;cn×1＝[c1c2…cn]T;ebn×1＝[eb1eb2…ebn]T;ecn×1＝[ec1ec2…ecn]T。

假定根據(jù)先驗信息,要求模型的截距參數(shù)β1和斜率參數(shù)β2必須在如下數(shù)值范圍內(nèi)

式(7)相當(dāng)于對平面擬合模型附加了4個參數(shù)β的不等式約束方程,相應(yīng)的G和z見表3。筆者共模擬了10個擬合點數(shù)據(jù),系數(shù)矩陣隨機量和觀測向量的中誤差見對應(yīng)觀測數(shù)據(jù)括號內(nèi)數(shù)值。估計結(jié)果見表4,可以看到,ICWTLS算法參數(shù)結(jié)果均滿足所有不等式約束條件(式(7)),而TLS解并不滿足不等式約束方程。即當(dāng)平差模型存在可靠的先驗信息時,將其作為附加約束條件,可得到滿足設(shè)定條件的估計結(jié)果。

表3 附有不等式約束的平面擬合模型數(shù)據(jù)Tab.3 Data set of the inequality-constrained plane fitting model

表4 附有不等式約束的平面擬合模型的ICWTLS估計結(jié)果Tab.4 Estimation of inequality-constrained WTLS and TLS of the adjustment model

5 結(jié) 論

當(dāng)EIV模型存在參數(shù)的先驗信息時,如模型某些參數(shù)取值應(yīng)在一定范圍內(nèi)或者地形擬合的邊界條件等就構(gòu)成了EIV模型的不等式約束條件。若平差計算顧及到約束條件,可以充分利用模型的先驗信息改善估計結(jié)果或者使得估計結(jié)果滿足設(shè)定的條件。目前僅有文獻(xiàn)[16—17]對ICTLS進(jìn)行了研究,但提出的算法只適用于系數(shù)矩陣全部元素隨機、觀測數(shù)據(jù)等權(quán)的特殊情況。此外,當(dāng)約束方程數(shù)量較大時,建立在窮舉法基礎(chǔ)上的ICTLS算法的計算量隨約束方程個數(shù)呈指數(shù)增長;當(dāng)約束方程數(shù)量過大時,算法甚至可能無法在有效時間內(nèi)進(jìn)行計算。本文將ICEIV模型改寫為更為一般化的ICPEIV模型,并且在整體最小二乘準(zhǔn)則下,將模型的求解轉(zhuǎn)化為標(biāo)準(zhǔn)的附有不等式約束的最優(yōu)化問題,提出的ICWTLS新算法能采用統(tǒng)一的方法估計結(jié)構(gòu)性系數(shù)矩陣、隨機和非隨機元素并存的各種情況,并且沒有限制觀測數(shù)據(jù)的精度和相關(guān)性。同時,算法的計算量不受約束方程個數(shù)的制約。論文通過兩個實例對ICWTLS新算法進(jìn)行了驗證和說明,計算結(jié)果表明算法較好地解決了現(xiàn)有算法的限制問題,新算法在實際應(yīng)用中簡單、有效、具有普遍適用性。

[1] XU P L,LIU J N,SHI C.Total Least Squares Adjustment inPartial Errors-in-variables Models:Algorithm and Statistical Analysis[J].Journal of Geodesy,2012,86(8):661-675.

[2] ADCOCK R J.Note on the Method of Least Squares[J].The Analyst,1877,4(6):183-184.

[3] PEARSON K.On Lines and Planes of Closest Fit to Systems of Points in Space[J].Philosophical Magazine Series 6, 1901,2(11):559-572.

[4] DEMING W E.The Application of Least Squares[J].Philosophical Magazine Series 7,1931,11(68):146-158.

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

[6] WASTON G A.Robust Counterparts of Errors-in-variables Problems[J].Computational Statistics&Data Analysis, 2007,52(2):1080-1089.

[7] VAN HUFFEL S,VANDEWALLE J.The Total Least Squares Problem:Computational Aspects and Analysis [M].Philadelphia:Society for Industrial and Applied Mathematics,1991.

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

[9] SCHAFFRIN B,LEE I,CHOI Y,et al.Total Least-squares (TLS)for Geodetic Straight-line and Plane Adjustment[J].Bollettino di Geodesia e Scienze Affini,2006,65(3):141-168.

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

[11] FANG X.Weighted Total Least Squares:Necessary and Sufficient Conditions,Fixed and Random Parameters[J].Journal of Geodesy,2013,87(8):733-749.

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

[13] SCHAFFRIN B,FELUS Y A.On Total Least-squares Adjustment with Constraints[C]∥International Association of Geodesy Symposia.Sapporo:International Association of Geodesy,2005,128:417-421.

[14] MAHBOUB V,SHARIFI M A.On Weighted Total Least Squares with Linear and Quadratic Constraints[J].Journal of Geodesy,2013,87(3):607-608.

[15] FANG Xing.A Structured and Constrained Total Least-Squares Solution with Cross-covariances[J].Studia Geophysica et Geodaetica,2014,58(1):1-16.

[16] DE MOOR B.Total Linear Least Squares with Inequality Constraints[R].Delft:Department of Electrical Engineering, 1990.

[17] ZH ANG Songlin,TONG Xiaohua,ZH ANG Kun.A Solution to EIV Model with Inequality Constraints and Its Geodetic Applications[J].Journal of Geodesy,2013, 87(1):89-99.

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

[19] NOCEDAL J,WRIGHT S J.Numerical Optimization[M].New York:Springer,2006.

[20] FLETCHER R.Practical Methods of Optimization[M].New York:John Wiley&Sons,2000.

(責(zé)任編輯:陳品馨)

Weighted Total Least Squares Algorithm with Inequality Constraints

ZENG Wenxian1,FANG Xing1,LIU Jingnan1,2,YAO Yibin1
1.School of Geodesy and Geomatics,Wuhan University,Wuhan 430079,China;2.Research Center of GNSS, Wuhan University,Wuhan 430079,China

Since the inequality-constrained total least squares(ICTLS)is strongly limited due to the combinational difficulty,adjustment of the partial errors-in-variables(EIV)model which is equipped with inequality constraints is investigated.The original ICTLS problem to a standard optimization problem is reconfigured in this paper,which can be solved by existing methods such as penalty based methods.The novel ICWTLS(inequality-constrained weighted TLS)algorithm can deal with the ICTLS problem with a structure coefficient matrix and a general weight matrix,and successfully avoid the combinatorial difficulty.The examples illustrate that the new algorithm proposed in this paper is efficient and simple, which can be used in a general case in practice.

total least squares;errors-in-variables model;inequality constraints;nonlinear program

ZENG Wenxian(1975—),female,PhD, majors in the theory and method of surveying data processing.

FANG Xing

P207

A

1001-1595(2014)10-1013-06

國家自然科學(xué)基金(41474006;41404005;41231174;41174012);中央高校基本科研基金(2042014kf053)

2013-12-17

曾文憲(1975—),女,博士,主要從事測量數(shù)據(jù)處理理論與應(yīng)用的研究。

E-mail:wxzeng＠sgg.whu.edu.cn

方興

E-mail:xfang＠sgg.whu.edu.cn

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):1013-1018.(曾文憲,方興,劉經(jīng)南,等.附有不等式約束的加權(quán)整體最小二乘算法[J].測繪學(xué)報,2014,43(10):1013-1018.)

10.13485/j.cnki.11-2089.2014.0173

修回日期:2014-08-29