張重陽，胡 川，李成洪，范小猛

(重慶交通大學(xué) 土木工程學(xué)院，重慶 400074)

整體最小二乘(total least squares，TLS)方法可以有效解決系數(shù)矩陣含有誤差的參數(shù)估計(jì)問題，得到了學(xué)者的廣泛關(guān)注。1980年，Golub和Van Loan從數(shù)值分析的角度用奇異值分解(singular value decomposition,SVD)方法得到了變量誤差(error-in-variable，EIV)模型的TLS解[1]，自此以后TLS法得到快速發(fā)展及廣泛應(yīng)用[2]。在測(cè)繪領(lǐng)域，通常會(huì)將TLS視作一個(gè)特殊的非線性問題，更多地是采用迭代解法進(jìn)行TLS平差。Schaffrin[3-4]為解決異方差和結(jié)構(gòu)EIV模型的參數(shù)估計(jì)問題，提出了加權(quán)整體最小二乘法(weighted total least squares，WTLS)。EIV模型在不同平差問題下具有不同的結(jié)構(gòu)特征[5]，部分變量誤差(partial error in variable，PEIV)模型[6-7,8]能有效解決不同結(jié)構(gòu)特征的EIV 模型參數(shù)估計(jì)問題。Fang[9]詳細(xì)地研究了WTLS解存在的充要條件和隨機(jī)參數(shù)估計(jì)問題，并闡述了不同迭代方法的異同。除此之外，也有研究人員從經(jīng)典平差的角度出發(fā)，提出采用標(biāo)準(zhǔn)最小二乘法(standard least squares，SLS)來求解EIV模型，得到WTLS解算結(jié)果(即SLS-WTLS)[10]，該法可直接利用間接平差的結(jié)果，有助于建立與經(jīng)典平差的聯(lián)系。WTLS迭代算法通常以加權(quán)最小二乘(weighted least squares，WLS)平差結(jié)果作為迭代初值，獲得局部最優(yōu)解，但當(dāng)初值離WLS解較遠(yuǎn)時(shí)，迭代解法就存在法矩陣奇異和算法發(fā)散的問題。

為了解決此問題，Schaffrin等[3]在進(jìn)行算法設(shè)計(jì)時(shí)，會(huì)預(yù)先計(jì)算一個(gè)近似WTLS解作為迭代初值，以確保算法具有良好的穩(wěn)定性及收斂性。另一方面，同倫方法也能解決非線性模型的初值依賴問題[11]。連續(xù)同倫法[12]起源于延拓法，主要用于求解方程組[13-15],但早期尚未形成穩(wěn)定有效算法。直到Li和Yorke提出預(yù)估-校正算法(Li-Yorke法)[16]，讓連續(xù)同倫法實(shí)現(xiàn)大范圍局部收斂。同倫方法在大地測(cè)量領(lǐng)域已有一些相關(guān)研究，如陶本藻等[17]將同倫方法引入GPS偽距定位模型的LS平差計(jì)算之中，提高估計(jì)精度；張勤等[18]在理論上將同倫方法與WLS相結(jié)合，解決了WLS法的秩虧平差問題；游為等[19]將非線性同倫最小二乘模型應(yīng)用于任意三維坐標(biāo)轉(zhuǎn)換，在初值離真值較遠(yuǎn)時(shí)獲得了七參數(shù)的穩(wěn)定解。前述研究的重點(diǎn)均在于非線性最小二乘(nonlinear least squares，NSL)，為了解決初值對(duì)SLS-WTLS法收斂性影響的問題，進(jìn)一步拓展EIV平差理論，文中以SLS-WTLS理論為媒介，采用同倫方法求解EIV模型，并提出同倫WTLS平差算法。

1 同倫WTLS

1.1 連續(xù)同倫方法及其解算

同倫來源于拓?fù)鋽?shù)學(xué)中的延拓概念，連續(xù)同倫方法是同倫方法中的一類，其算法是由Chow等提出的[15]。假設(shè)要求解的模型為y=F(x)，為方程組yi=f(xi),(i=1,…,2，n)對(duì)應(yīng)的向量值函數(shù)，同倫函數(shù)[19]可以定義為h(t,x)∈Rn+1：

(1)

其中

(2)

(3)

式中:B,G分別是F(X)的一階(Jacobian matrix)和二階偏導(dǎo)(Hessian matrix)矩陣?？梢圆捎妙A(yù)估-校正法對(duì)式(3)進(jìn)行解算，該算法具體實(shí)現(xiàn)過程可以參考文獻(xiàn)[16]。

1.2 同倫WTLS平差模型及其解算

經(jīng)研究，當(dāng)以WLS解作為初值時(shí)，SLS-WTLS法能快速收斂并取得理想估計(jì)結(jié)果，但在初值精度較差的情況下，該算法會(huì)存在奇異和發(fā)散問題。為此文中引入同倫思想對(duì)該問題予以解決，一般的EIV函數(shù)模型可以表達(dá)為：

y+ey=(A+EA)ξ.

(4)

其中，y=(y1,y2,…,yn)T為觀測(cè)向量，ξ=(ξ1,…,ξm)T為參數(shù)向量，A為n×m階系數(shù)矩陣，ey與EA分別代表y與A中的偶然誤差。

EIV隨機(jī)模型表達(dá)為：

(5)

QA=Q0?Qx.

(6)

(7)

隨機(jī)模型為：

(8)

由法方程得到的參數(shù)計(jì)算式為：

(9)

根據(jù)前述連續(xù)同倫方法對(duì)式(7)進(jìn)行解算，首先按照標(biāo)準(zhǔn)最小二乘準(zhǔn)則得到正交性條件方程為：

(10)

(11)

(12)

對(duì)式(12)求全微分，得到WTLS模型為：

(13)

進(jìn)一步可以得到，

(14)

由式(13)計(jì)算切向量vk，即：

(15)

預(yù)估步長常用計(jì)算方法有歐拉預(yù)估法及四階龍格-庫塔(Runge-Kutta)法等[21]，文中采用歐拉預(yù)估法。由式(16)可進(jìn)一步得到預(yù)估值(tk+1,ξk+1)，即：

(tk+1,ξk+1)=(tk,ξk)+Δt·vk.

(16)

(H(tk+1,ξk+1)).

(17)

根據(jù)預(yù)估-校正法的原理，對(duì)同倫解曲線λ(s)=(t(s),ξ(s))進(jìn)行跟蹤，從t=0和ξ=ξ0端點(diǎn)開始，直到t=1且g2(ξ*)=0，解ξ*即為求得的參數(shù)估計(jì)值。外符合精度加權(quán)殘差平方和(TSSR)及內(nèi)符合精度驗(yàn)后單位權(quán)方差計(jì)算式為：

(18)

(19)

根據(jù)SLS理論，參數(shù)估計(jì)值的協(xié)因數(shù)陣為：

(20)

(21)

2 同倫WTLS平差算法

根據(jù)同倫方法的基本原理，推導(dǎo)得到同倫WTLS平差計(jì)算式，為了得到穩(wěn)定有效的結(jié)果，算法流程設(shè)計(jì)如下：

輸入觀測(cè)向量，系數(shù)矩陣A及對(duì)應(yīng)的權(quán)陣Px，Py。

1)初始化：確定初始點(diǎn)(t(0),ξ(0))=(0,ξLS)，并給定步長Δt與閾值ε1=10-7、ε2=10-5，ε3取值與步長有關(guān),k=0。

3)校正：按式(17)以(tk+1,ξk+1)為初始值，采用牛頓迭代法得到tk+1參數(shù)校正序列{ξk+1(i)},i=1,2,…。為了與算法計(jì)算次數(shù)k有所區(qū)分，以i值代表牛頓迭代計(jì)算次數(shù)。牛頓迭代的終止條件[22]為:

(22)

4)檢驗(yàn)：判斷tk+1，若|tk+1-1|<ε3，輸出ξk+1*并以式(18)～式(21)進(jìn)行精度評(píng)定，否則以(tk+1,ξk+1(i))為初始值，轉(zhuǎn)至步驟(2)。

3 實(shí)驗(yàn)分析

為了驗(yàn)證所提算法的可行性，評(píng)估同倫加權(quán)整體最小二乘算法的性能，本節(jié)以直線擬合和平面坐標(biāo)變換問題為例進(jìn)行討論。

3.1 直線擬合

直線坐標(biāo)數(shù)據(jù)來源于文獻(xiàn)[23]，對(duì)應(yīng)的權(quán)值來源于文獻(xiàn)[24]。為了方便比較，文中將文獻(xiàn)[9]的算法簡(jiǎn)稱為Fang-WTLS，并將其與同倫WTLS法分別對(duì)前述直線數(shù)據(jù)進(jìn)行擬合，3種方法的參數(shù)估計(jì)結(jié)果見表1。

表1 同倫加權(quán)整體最小二乘算法參數(shù)估計(jì)結(jié)果

經(jīng)研究計(jì)算發(fā)現(xiàn)，文獻(xiàn)[3,10,25]所提算法的計(jì)算結(jié)果與Fang-WTLS法相一致。根據(jù)表1可知，在步長選取恰當(dāng)?shù)臈l件下，新算法能夠得到可靠穩(wěn)定的結(jié)果，與Fang-WTLS法保持一致，由此可見，同倫WTLS法具有可行性。另一方面，隨著步長的減小，同倫解曲線的零點(diǎn)數(shù)量在逐漸增加，其收斂速度亦隨之減小。

在此基礎(chǔ)之上，針對(duì)SLS-WTLS算法自身存奇異及發(fā)散問題進(jìn)一步研究。為此通過不斷增大新初值與WLS解的距離‖Δξ0‖=‖ξ0-ξWLS‖，對(duì)比兩種算法的收斂性，表2列出了兩種算法收斂性差異。

表2 兩種算法收斂性比較

當(dāng)‖Δξ0‖逐漸增大時(shí)，SLS-WTLS法在進(jìn)行迭代的過程中出現(xiàn)了法方程奇異問題和發(fā)散問題，而同倫WTLS算法在步長取0.000 01時(shí)估計(jì)結(jié)果穩(wěn)定收斂并且精度可達(dá)mm級(jí)以上。

綜上所述，同倫WTLS法具有可行性，在步長選取恰當(dāng)?shù)臈l件下，得到的結(jié)果與Fang-WTLS法相一致。在初值離真值較遠(yuǎn)時(shí)，新算法仍能夠得到穩(wěn)定收斂解。

3.2 二維坐標(biāo)變換

二維坐標(biāo)變換問題的核心是確定目標(biāo)坐標(biāo)和原坐標(biāo)系之間變換的參數(shù)。算例的實(shí)驗(yàn)數(shù)據(jù)來源于文獻(xiàn)[26]，在步長分別取0.1、0.01和0.001時(shí)采用同倫WTLS法得到4項(xiàng)參數(shù)估計(jì)結(jié)果及精度列于表3。

根據(jù)表3可知，Δt取0.001時(shí)，平面坐標(biāo)變換的參數(shù)估計(jì)結(jié)果及精度指標(biāo)相對(duì)于WLS法來說，更符合實(shí)際的情況，再次證明了同倫WTLS算法的可行性。圖1描繪了各項(xiàng)參數(shù)估計(jì)結(jié)果隨步長取值變化情況。

從圖1可知，當(dāng)步長較大時(shí)，收斂速度較慢，t接近0.8左右時(shí)曲線才接近平緩，這是由于誤差累積造成的。而步長取0.001時(shí)解曲線在t接近0.1時(shí)就接近收斂狀態(tài)，最終結(jié)果與Fang-WTLS法的估計(jì)結(jié)果及精度指標(biāo)相一致。

為了進(jìn)一步驗(yàn)證新算法的收斂性及穩(wěn)定性，在初值取值不同的情況下，采用同倫WTLS與SLS-WTLS法參數(shù)估計(jì)結(jié)果及精度見表4。

根據(jù)表4可知，就坐標(biāo)轉(zhuǎn)換問題而言，當(dāng)初始距真值較遠(yuǎn)時(shí)，SLS-WTLS法在迭代時(shí)會(huì)存在矩陣計(jì)算奇異而不收斂，由此導(dǎo)致最終結(jié)果不準(zhǔn)確或是計(jì)算失敗。二維坐標(biāo)轉(zhuǎn)換算例結(jié)果表明同倫WTLS法具有良好的收斂性，這與前述直線擬合算例是一致的。

表3 二維坐標(biāo)變換問題的同倫WTLS法參數(shù)估計(jì)結(jié)果

表4 不同初值下兩種算法平面坐標(biāo)變換參數(shù)估計(jì)值

圖1 不同步長取值對(duì)應(yīng)坐標(biāo)變換各項(xiàng)參數(shù)估計(jì)結(jié)果

4 結(jié)束語

文中以SLS-WTLS迭代解法存在的初值問題為出發(fā)點(diǎn)，將同倫方法與加權(quán)整體最小二乘法相結(jié)合，提出一種用于線性EIV模型參數(shù)估計(jì)新方法。

1)在步長選取適當(dāng)?shù)臈l件下，同倫WTLS方法估計(jì)結(jié)果能與其他WTLS算法保持一致，新算法具有可行性。

2)同倫WTLS算法具有大范圍收斂的性質(zhì)，解決了SLS-WTLS法在近似值距真值較遠(yuǎn)時(shí)存在法矩陣奇異和算法發(fā)散問題。

3)參數(shù)估計(jì)及精度與步長取值有著較為密切的關(guān)系，選取恰當(dāng)?shù)牟介L，參數(shù)估計(jì)精度就越高，結(jié)果也就越可靠。