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

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

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

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

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

1 非線性定位方程及其解

非線性測(cè)距定位方程可以表示為

L=d(x)+ε

(1)

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

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

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

(2)

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

dL=Jdx+δ

(3)

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

JTP(JTdx-dL)=0

(4)

即

(5)

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

(6)

(7)

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

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

(8)

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

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

(9)

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

梯度法收斂速度較慢，但優(yōu)點(diǎn)是不需要矩陣求逆運(yùn)算。目前最常用的迭代算法有牛頓法、高斯-牛頓法和各種改進(jìn)的牛頓法等。其中，高斯-牛頓法的基本思想就是在初值x0處對(duì)非線性模型進(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)度大，測(cè)距方程出現(xiàn)病態(tài)問(wèn)題時(shí)，高斯-牛頓法可能會(huì)出現(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)造一個(gè)迭代序列，使得這個(gè)解序列滿足收斂條件。但是不同算法的收斂性存在差別，尤其是局部收斂性。當(dāng)方程呈現(xiàn)病態(tài)性時(shí)，會(huì)直接影響迭代算法的穩(wěn)定性，對(duì)解的收斂性產(chǎn)生一定的影響。本文將重點(diǎn)討論方程的多解性問(wèn)題。

2 非線性定位方程的多解性問(wèn)題

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

2.1 多解性及辨識(shí)

對(duì)于測(cè)距方程，當(dāng)n=3時(shí)，由文獻(xiàn)[15]可知，肯定存在兩個(gè)解。共線時(shí)，存在無(wú)窮解。可采用解析法進(jìn)行分析。

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

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

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

3 算 例

3.1 算例1

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

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

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

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

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

3.2 算例2

當(dāng)n≥4時(shí)，選取P1、P2、…、P9中所有的已知點(diǎn)及這9個(gè)已知點(diǎn)到未知點(diǎn)P10的觀測(cè)距離，算例中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)了波動(dòng)，最后穩(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é)果與真值很近似，說(shuō)明在X方向上收斂性較好，因而下面討論初值取不同值時(shí)，固定初值的X大小(本文設(shè)X=0)，Y和Z取值范圍為[-10:1:10],收斂結(jié)果見(jiàn)表3。

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

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

圖3 空間位置初值對(duì)收斂解的影響

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

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

4 結(jié) 語(yǔ)

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

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

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