曾文憲,方 興,劉經(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

針對(duì)現(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)陣限制條件擴(kuò)展到了一般性權(quán)矩陣,將要求系數(shù)矩陣元素全部隨機(jī)的限定條件擴(kuò)展到了可同時(shí)包含隨機(jī)和非隨機(jī)元素的一般情況,并且新算法解決了現(xiàn)有算法計(jì)算量受制于約束方程數(shù)量的缺陷。實(shí)例計(jì)算表明,本文提出的算法簡單、有效,具有普遍適用性。

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

1 引 言

整體最小二乘估計(jì)(total least squares, TLS)作為EIV[1](errors-in-variables)模型的嚴(yán)密估計(jì)方法,目前已廣泛應(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)計(jì)意義上的TLS算法。文獻(xiàn)[5]提出了著名的奇異值分解(singular value decomposition,SVD)算法。文獻(xiàn)[6]研究了穩(wěn)健整體最小二乘算法(robust TLS)等。文獻(xiàn)[7]證明了當(dāng)觀測數(shù)趨于無窮大時(shí),TLS估計(jì)具有漸進(jìn)無偏性。大地測量領(lǐng)域針對(duì)普遍存在的觀測值不等精度、相關(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ù)估計(jì)存在先驗(yàn)信息時(shí),EIV模型擴(kuò)展為附有約束的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)。

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

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

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

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

式中,y表示n×1的觀測向量;ey表示y的n×1隨機(jī)誤差向量,且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隨機(jī)誤差矩陣, 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中元素全部為隨機(jī)量。當(dāng)A中存在非隨機(jī)的固定元素或者呈現(xiàn)結(jié)構(gòu)性特征時(shí),必須對(duì)系數(shù)矩陣的協(xié)因數(shù)陣[18]或者對(duì)系數(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)勢(shì)主要體現(xiàn)為:①將ICEIV模型中系數(shù)矩陣元素全部為隨機(jī)量的限定擴(kuò)展到了同時(shí)包含隨機(jī)和非隨機(jī)元素的一般情況;②ICPEIV模型可解析結(jié)構(gòu)性系數(shù)矩陣,保證了算法的統(tǒng)一性;③ICPEIV模型A中參與平差計(jì)算的待估量個(gè)數(shù)m小于等于ICEIV中A中元素個(gè)數(shù)nt,尤其當(dāng)系數(shù)矩陣中的隨機(jī)量個(gè)數(shù)較少時(shí),ICPEIV模型形式可以大大減少待估計(jì)量。因此,ICPEIV模型更具一般性,以下在不限定觀測數(shù)據(jù)權(quán)陣的等精度和相關(guān)性的一般情況下,即ICPEIV模型(2)中權(quán)矩陣ω和w為任意正定對(duì)稱陣,筆者提出了普遍適用的附有不等式約束的加權(quán)整體最小二乘新算法。

3 (ICWTLS)算法

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

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

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

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

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

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

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

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

4 實(shí)例分析

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

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

實(shí)例1的模擬數(shù)據(jù)見表1,模型包含7×1待估參數(shù)向量β以及70×1待估系數(shù)向量ˉa,系數(shù)矩陣A的全部元素為隨機(jī)量,A和y中所有觀測元素不相關(guān)且中誤差均為0.1,其中參數(shù)真值β～見表2。為了改進(jìn)模型的估計(jì)結(jié)果,利用參數(shù)的先驗(yàn)信息設(shè)計(jì)了18個(gè)不等式約束方程(見表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)從算法的計(jì)算量比較,本文算法只需迭代計(jì)算23次,而窮舉法理論上排列組合所需計(jì)算次數(shù)為218≈260 000。當(dāng)模型約束方程個(gè)數(shù)較多的情況下,基于窮舉法的ICTLS算法計(jì)算量要遠(yuǎn)大于本文提出的ICWTLS算法。當(dāng)不等式約束方程過多時(shí),基于窮舉法的算法甚至無法在有效時(shí)間內(nèi)求解。

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

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

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

由上式可以看到平面擬合模型的系數(shù)矩陣同時(shí)包含了非隨機(jī)和隨機(jī)元素,即第1列元素為已知量,第2列和第3列是由觀測數(shù)據(jù)構(gòu)成的隨機(jī)量,現(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ù)先驗(yàn)信息,要求模型的截距參數(shù)β1和斜率參數(shù)β2必須在如下數(shù)值范圍內(nèi)

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

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

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

5 結(jié) 論

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

10.13485/j.cnki.11-2089.2014.0173

修回日期:2014-08-29