同倫加權整體最小二乘平差算法

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

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

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

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

1 同倫WTLS

1.1 連續同倫方法及其解算

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

(1)

其中

(2)

(3)

式中:B,G分別是F(X)的一階(Jacobian matrix)和二階偏導(Hessian matrix)矩陣。可以采用預估-校正法對式(3)進行解算，該算法具體實現過程可以參考文獻[16]。

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

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

y+ey=(A+EA)ξ.

(4)

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

EIV隨機模型表達為：

(5)

QA=Q0?Qx.

(6)

(7)

隨機模型為：

(8)

由法方程得到的參數計算式為：

(9)

根據前述連續同倫方法對式(7)進行解算，首先按照標準最小二乘準則得到正交性條件方程為：

(10)

(11)

(12)

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

(13)

進一步可以得到，

(14)

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

(15)

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

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

(16)

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

(17)

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

(18)

(19)

根據SLS理論，參數估計值的協因數陣為：

(20)

(21)

2 同倫WTLS平差算法

根據同倫方法的基本原理，推導得到同倫WTLS平差計算式，為了得到穩定有效的結果，算法流程設計如下：

輸入觀測向量，系數矩陣A及對應的權陣Px，Py。

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

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

(22)

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

3 實驗分析

為了驗證所提算法的可行性，評估同倫加權整體最小二乘算法的性能，本節以直線擬合和平面坐標變換問題為例進行討論。

3.1 直線擬合

直線坐標數據來源于文獻[23]，對應的權值來源于文獻[24]。為了方便比較，文中將文獻[9]的算法簡稱為Fang-WTLS，并將其與同倫WTLS法分別對前述直線數據進行擬合，3種方法的參數估計結果見表1。

表1 同倫加權整體最小二乘算法參數估計結果

經研究計算發現，文獻[3,10,25]所提算法的計算結果與Fang-WTLS法相一致。根據表1可知，在步長選取恰當的條件下，新算法能夠得到可靠穩定的結果，與Fang-WTLS法保持一致，由此可見，同倫WTLS法具有可行性。另一方面，隨著步長的減小，同倫解曲線的零點數量在逐漸增加，其收斂速度亦隨之減小。

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

表2 兩種算法收斂性比較

當‖Δξ0‖逐漸增大時，SLS-WTLS法在進行迭代的過程中出現了法方程奇異問題和發散問題，而同倫WTLS算法在步長取0.000 01時估計結果穩定收斂并且精度可達mm級以上。

綜上所述，同倫WTLS法具有可行性，在步長選取恰當的條件下，得到的結果與Fang-WTLS法相一致。在初值離真值較遠時，新算法仍能夠得到穩定收斂解。

3.2 二維坐標變換

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

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

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

為了進一步驗證新算法的收斂性及穩定性，在初值取值不同的情況下，采用同倫WTLS與SLS-WTLS法參數估計結果及精度見表4。

根據表4可知，就坐標轉換問題而言，當初始距真值較遠時，SLS-WTLS法在迭代時會存在矩陣計算奇異而不收斂，由此導致最終結果不準確或是計算失敗。二維坐標轉換算例結果表明同倫WTLS法具有良好的收斂性，這與前述直線擬合算例是一致的。

表3 二維坐標變換問題的同倫WTLS法參數估計結果

表4 不同初值下兩種算法平面坐標變換參數估計值

圖1 不同步長取值對應坐標變換各項參數估計結果

4 結束語

文中以SLS-WTLS迭代解法存在的初值問題為出發點，將同倫方法與加權整體最小二乘法相結合，提出一種用于線性EIV模型參數估計新方法。

1)在步長選取適當的條件下，同倫WTLS方法估計結果能與其他WTLS算法保持一致，新算法具有可行性。

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

3)參數估計及精度與步長取值有著較為密切的關系，選取恰當的步長，參數估計精度就越高，結果也就越可靠。