測距定位方程的多解性及其非線性最小二乘迭代算法

齊 珂，曲國慶，薛樹強(qiáng)，劉以旭，楊文龍，韓德強(qiáng)

(1. 山東理工大學(xué)，山東 淄博 255049； 2. 中國測繪科學(xué)研究院，北京 100830)

大地測量觀測模型一般為非線性模型[1]，運(yùn)用較廣的經(jīng)典平差中均為對非線性觀測方程在參數(shù)的近似值處進(jìn)行泰勒級數(shù)展開，取一階項(xiàng)將其轉(zhuǎn)化為線性方程進(jìn)行平差計(jì)算[2]。但是這要求未知數(shù)近似值與真值充分接近且非線性模型的非線性程度較弱才能進(jìn)行線性化[3]。固有曲率和參數(shù)效應(yīng)是刻畫非線性模型線性化的數(shù)量指標(biāo)，能夠評估統(tǒng)計(jì)推斷效果的優(yōu)劣程度、適用條件和容許誤差[4-6]。

當(dāng)測距方程的非線性強(qiáng)度較高時，傳統(tǒng)的平差方法將不再適用[7]。在非線性模型參數(shù)估計(jì)方面，高斯-牛頓法、最速下降法、阻尼最小二乘法[8-9]等研究較多。文獻(xiàn)[10]針對非線性整體最小二乘問題給出了相應(yīng)的迭代算法，推導(dǎo)了基于牛頓法的平差迭代新公式。在參數(shù)估計(jì)的直接解法方面，文獻(xiàn)[11]給出了顧及二次項(xiàng)和三次項(xiàng)影響的直接解。然而這些研究只是針對如何求解這一問題進(jìn)行研究，對求解的多樣性問題研究較少。

本文在非線性測距方程的迭代算法和直接解法的基礎(chǔ)上，針對非線性方程的多解性進(jìn)行討論，并且根據(jù)解的收斂性對目前常用的迭代算法進(jìn)行分析，討論不同迭代算法的特點(diǎn)，最后證實(shí)封閉牛頓迭代算法在方程存在多解時具有更好的局部收斂性。

1 非線性定位方程及其解

非線性測距定位方程可以表示為

L=d(x)+ε

(1)

式中，L為觀測向量(L1,L2,…,Ln)；d(x)為未知點(diǎn)到已知點(diǎn)的歐氏距離；x為未知點(diǎn)坐標(biāo)；ε為觀測誤差(ε1,ε2,…,εn)。

式(1)為非線性方程，將其在初值x0處線性化并代入式(1)中，可得線性化方程組

L=Jdx+r+ε+d(x0)

(2)

式中，r為線性化殘余項(xiàng)。上式可改寫為

dL=Jdx+δ

(3)

式中，J為雅克比矩陣；δ=r+ε。

JTP(JTdx-dL)=0

(4)

即

(5)

此時，可以得出未知點(diǎn)的坐標(biāo)為

(6)

(7)

則可得以下非線性正交條件

h(x)=J(x)TPV(x)=0

(8)

式(8)為一相容方程組，此相容方程的解即為定位方程的非線性最小二乘解。

在求解非線性方程(8)時，一般需要迭代算法進(jìn)行迭代計(jì)算。最簡單的迭代方法是梯度法，也稱最速下降法，它是通過沿梯度下降的方向求解極小值的方法確定方程的解，以下迭代公式可確保算法收斂[13]

(9)

式中，tr( )為矩陣的跡；P為觀測權(quán)陣。

梯度法收斂速度較慢，但優(yōu)點(diǎn)是不需要矩陣求逆運(yùn)算。目前最常用的迭代算法有牛頓法、高斯-牛頓法和各種改進(jìn)的牛頓法等。其中，高斯-牛頓法的基本思想就是在初值x0處對非線性模型進(jìn)行線性近似，按傳統(tǒng)的平差方法求出近似值，然后反復(fù)迭代，直至滿足迭代條件，其迭代公式為

xk+1=xk+N-1(xk)h(xk)

(10)

式中，N(x)=J(x)TPJ(x)，高斯-牛頓法的適用條件為初值精度足夠高、殘差較小、非線性強(qiáng)度較弱。當(dāng)初值遠(yuǎn)離真值或非線性強(qiáng)度大，測距方程出現(xiàn)病態(tài)問題時，高斯-牛頓法可能會出現(xiàn)不收斂。

本文采用的另一種迭代算法是一種改進(jìn)的牛頓算法，為封閉牛頓法。由文獻(xiàn)[14]可知，求解方程(4)的封閉牛頓法的迭代計(jì)算公式為

xk+1=xk+H-1(xk)h(xk)

(11)

式中，H(x)為x的海森矩陣，表示為

(12)

需要進(jìn)一步指出，求解非線性方程的迭代算法不同，但其最終收斂的解一般卻相同。迭代算法的理論基礎(chǔ)為壓縮影射原理，即構(gòu)造一個迭代序列，使得這個解序列滿足收斂條件。但是不同算法的收斂性存在差別，尤其是局部收斂性。當(dāng)方程呈現(xiàn)病態(tài)性時，會直接影響迭代算法的穩(wěn)定性，對解的收斂性產(chǎn)生一定的影響。本文將重點(diǎn)討論方程的多解性問題。

2 非線性定位方程的多解性問題

非線性方程組h(x)=JTPV=0的解可能不唯一，并且不同的解對應(yīng)的VTPV極值點(diǎn)也是存在差異的。從統(tǒng)計(jì)意義上講，也并非VTPV數(shù)值最小的解對應(yīng)最佳參數(shù)估計(jì)。在實(shí)際中，希望與各類先驗(yàn)信息吻合的解，包括先驗(yàn)初值、先驗(yàn)觀測方差等。

2.1 多解性及辨識

對于測距方程，當(dāng)n=3時，由文獻(xiàn)[15]可知，肯定存在兩個解。共線時，存在無窮解?？刹捎媒馕龇ㄟM(jìn)行分析。

當(dāng)n≥4時，解的情況將變得復(fù)雜，本文將采用常用的幾種迭代算法對方程進(jìn)行迭代求解并計(jì)算出相應(yīng)的VTPV的數(shù)值，結(jié)合先驗(yàn)信息對收斂解及解的VTPV數(shù)值進(jìn)行分析。筆者建議，通過后驗(yàn)單位權(quán)方差與實(shí)際觀測數(shù)據(jù)的精度信息進(jìn)行對照，得到具有實(shí)際意義的解。

2.2 附加穩(wěn)定泛函的正則化算法

測距方程出現(xiàn)病態(tài)問題時，最常用的方法除了解析法，還有直接法。為了改善雅克比矩陣秩虧或病態(tài)對參數(shù)估值的影響，一般情況為利用先驗(yàn)信息對其附加強(qiáng)約束，即引入穩(wěn)定泛函極小準(zhǔn)則。

3 算 例

3.1 算例1

本文采用文獻(xiàn)[16]模擬的測邊網(wǎng)，P1、P2、…、P9為9個已知點(diǎn)，以及這9個已知點(diǎn)到未知點(diǎn)P10(模擬值為(0,0,0))的觀測距離(如圖1所示)，控制點(diǎn)坐標(biāo)和觀測值見表1，距離為等精度觀測，中誤差為0.001 m，根據(jù)9個觀測距離確定未知點(diǎn)的坐標(biāo)。

圖1 病態(tài)測邊網(wǎng)

點(diǎn)號坐標(biāo)XYZ觀測距離P123.00010.0000.01025.0786P2-10.0009.9900.00014.1345P335.00010.010-0.01036.4159P4100.00019.9900.005101.4794P5-36.00010.005-0.00037.3642P60.00010.010-0.00510.0100P756.0009.9950.01056.9961P8-15.00010.015-0.01018.0359P9-1.700010.0080.01510.1506

當(dāng)n=3時，選取P2、P3、P4這3個已知點(diǎn)及對應(yīng)的觀測距離，算例中N=JTJ的條件數(shù)為4.585 1×103，嚴(yán)重病態(tài)。計(jì)算的坐標(biāo)近似值為(0，0，0)，分別采用解析法和直接法進(jìn)行定位解算(直接法采用Matlab中solve函數(shù)解非線性方程組)，結(jié)果見表2。

從上述3種方法的計(jì)算結(jié)果可以看出，定位方程至少是存在3個解的。不同算法的解是不同的，根據(jù)先驗(yàn)信息(未知點(diǎn)模擬值為(0,0,0))可以看出，只有封閉牛頓法第3個解比較符合，但是其VTPV數(shù)值是最大的,這進(jìn)一步說明非線性最小二乘求解是求解區(qū)域VTPV數(shù)值最小，而不是全局最小。根據(jù)VTPV數(shù)值的大小判斷解的情況似乎不是很合理，目前來講只有通過先驗(yàn)信息來判斷解是否可行。但是，可以看出，當(dāng)其他方法與封閉牛頓法不同時，封閉牛頓法更接近真值。

表2 不同算法的數(shù)值收斂解

3.2 算例2

當(dāng)n≥4時，選取P1、P2、…、P9中所有的已知點(diǎn)及這9個已知點(diǎn)到未知點(diǎn)P10的觀測距離，算例中N=JTJ的條件數(shù)為4.585 1×103，嚴(yán)重病態(tài)。計(jì)算的坐標(biāo)近似值為(0，0，0)。分別采用封閉牛頓法和高斯-牛頓法進(jìn)行計(jì)算未知點(diǎn)P10，得到的點(diǎn)位坐標(biāo)迭代解序列如圖2所示。

圖2 點(diǎn)位坐標(biāo)迭代解序列

圖2給出了封閉牛頓法和高斯-牛頓法的點(diǎn)位坐標(biāo)迭代解序列。圖中，縱坐標(biāo)為收斂坐標(biāo)值，橫軸為迭代次數(shù)。從圖中可以看出，高斯-牛頓法收斂速度較慢，在收斂中出現(xiàn)了波動，最后穩(wěn)定收斂于(-0.027,5.366 3,8.877 8)，但與真值(0,0,0)相差很大；封閉牛頓法收斂速度較快，迭代解序列也很穩(wěn)定，并且最終收斂于(0.063，0.025，0.017)，與真值(0,0,0)十分接近。可以看出，在收斂速度與穩(wěn)定性上，封閉牛頓法更加優(yōu)越一些。

可以看出，封閉牛頓法與高斯-牛頓法在X方向上收斂結(jié)果與真值很近似，說明在X方向上收斂性較好，因而下面討論初值取不同值時，固定初值的X大小(本文設(shè)X=0)，Y和Z取值范圍為[-10:1:10],收斂結(jié)果見表3。

表3 不同算法的迭代收斂解

從表中可以看出，給出任意初值，方程都能收斂到這幾個解，說明此非線性方程至少有3個解，封閉牛頓法與高斯-牛頓法收斂到相同的兩個解，并且封閉牛頓法會多收斂到一個更接近符合條件的解。其原因是：在局部收斂性上，封閉牛頓法比高斯-牛頓法要好，高斯-牛頓法無法收斂到解(0.063 3,0.025 5,0.017 7)。為能更進(jìn)一步說明這3個值為非線性方程的解，本文選取了其中一個值(0.063 3,0.025 5,0.017 7)附近的VTPV圖，如圖3所示。

圖3 空間位置初值對收斂解的影響

圖3中，圖中心點(diǎn)處是解(0.063 3,0.025 5,0.017 7)的VTPV值，顏色深淺代表取不同Y和Z時VTPV值的大小，從圖3中可以看出，中心點(diǎn)即收斂解處的VTPV的值是極小值，進(jìn)一步證實(shí)了這3個值即為非線性方程的解，同時也說明了封閉牛頓法是局部收斂的。

同時，在圖1的對比中可以看出，封閉牛頓法的迭代次數(shù)明顯比高斯牛頓法少，封閉牛頓法明顯更加有效率；從表1中可以看出，當(dāng)取不同的初值時，封閉牛頓法與高斯-牛頓法收斂到不同的解，因而，初值的選取對求解非線性方程很重要。同時也可以看出，兩種方法會收斂到不同的解，當(dāng)兩種方法的收斂解不同時，封閉牛頓法收斂到的解明顯更接近真值。這兩個試驗(yàn)結(jié)果都顯示了封閉牛頓法收斂到的VTPV數(shù)值較大的解為正確解，這是試驗(yàn)的偶然性還是有某種關(guān)系還需要進(jìn)一步試驗(yàn)分析。

4 結(jié) 語

當(dāng)非線性方程呈現(xiàn)病態(tài)性時，方程的非線性強(qiáng)度會很大，本文試驗(yàn)驗(yàn)證了解析法、高斯-牛頓法和封閉牛頓法的性能和收斂性質(zhì)，高斯-牛頓法明顯收斂速度較慢，并且會出現(xiàn)不穩(wěn)定的情況，試驗(yàn)表明封閉牛頓法的局部收斂性好于高斯-牛頓法。

當(dāng)非線性強(qiáng)度不高時，高斯-牛頓法與封閉牛頓法是沒有區(qū)別的，會收斂到相同的結(jié)果，解也是唯一的。當(dāng)非線性方程出現(xiàn)病態(tài)問題時會出現(xiàn)多解情況，解析法與高斯-牛頓法會得到兩個解，封閉牛頓法會得到3個解，并且其中兩個解與前兩種算法相同，但是這兩個解與真值相差甚遠(yuǎn)。封閉牛頓法多出的這個解與真值非常接近，在沒有先驗(yàn)信息的情況下如何區(qū)別這個解，目前還沒有方法。同時，這些解都為局部最優(yōu)解。

最后試驗(yàn)發(fā)現(xiàn)封閉牛頓法收斂到的正確解的VTPV數(shù)值是3個解中最大的，這可以為區(qū)別3個解中的正確解提供一個思路。