非線性總體最小二乘反演道路圓曲線參數(shù)算法

胡 川， 鄧明鏡， 趙立都

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

0 引 言

直線、緩和曲線和圓曲線是道路平面線形的基本元素。線形參數(shù)的反演，在公路普查、建立現(xiàn)有道路路線圖形數(shù)據(jù)庫等方面具有非常重要的作用，特別是圓曲線參數(shù)，影響著緩和曲線參數(shù)、道路里程的計(jì)算結(jié)果[1]。重心坐標(biāo)法[2-3]、附有參數(shù)的條件平差法[4]和最小二乘法[5]是3種較為常見的圓曲線參數(shù)估計(jì)方法。重心坐標(biāo)法是根據(jù)不在同一條直線上的3點(diǎn)確定一個(gè)圓，將所有觀測坐標(biāo)進(jìn)行排列組合，計(jì)算各組合點(diǎn)的圓心坐標(biāo)和半徑，取所有組合計(jì)算結(jié)果的平均值作為圓曲線參數(shù)的估計(jì)值。該方法的優(yōu)點(diǎn)是計(jì)算簡單，缺點(diǎn)是平差理論不夠嚴(yán)密；附有參數(shù)的條件平差是將圓曲線參數(shù)估計(jì)的非線性最小二乘模型轉(zhuǎn)化為附有參數(shù)的條件平差模型，采用經(jīng)典平差結(jié)果直接計(jì)算參數(shù)估計(jì)值。該方法的優(yōu)點(diǎn)是平差理論相對嚴(yán)密，缺點(diǎn)是參數(shù)估計(jì)結(jié)果易受參數(shù)近似值(初值)的影響[6]。最小二乘法以觀測點(diǎn)到圓曲線距離的平方和為極小值條件，采用換元法將其轉(zhuǎn)化為普通最小二乘平差模型[7]。該方法的優(yōu)點(diǎn)是避免計(jì)算非線性問題，缺點(diǎn)是其極小值條件不完全合理，且不具有明確的幾何意義。道路圓曲線參數(shù)反演既要計(jì)算參數(shù)[8-9]，還要估計(jì)誤差值，為剔除粗差影響提供基礎(chǔ)保障[10]。近年來，非線性總體最小二乘(nonlinear total least squares，NTLS)引起了測繪數(shù)據(jù)處理領(lǐng)域的研究興趣[11]，針對上述問題，筆者提出一種基于NTLS的道路圓曲線參數(shù)反演算法，直接用觀測誤差定義平差極小值條件(平差準(zhǔn)則)，減小參數(shù)初值對估計(jì)結(jié)果的影響，以提高道路圓曲線參數(shù)反演精度和理論的嚴(yán)密性。

1 NTLS反演道路圓曲線參數(shù)模型

用離散觀測數(shù)據(jù)反算道路線形參數(shù)，首先需要識別哪些點(diǎn)屬于哪類線性，例如直線、圓曲線；其次是選用合理、可靠的數(shù)學(xué)方法估計(jì)其最佳參數(shù)值。筆者假設(shè)已識別出圓曲線上的觀測數(shù)據(jù)，重點(diǎn)討論估計(jì)圓曲線參數(shù)的NTLS方法。

假設(shè)現(xiàn)有GNSS接收機(jī)測量的m個(gè)圓曲線路面段離散坐標(biāo) (xi，yi)，i=1,2,…,m。設(shè)該圓曲線段的圓心坐標(biāo)為(a，b)，半徑為r，則估計(jì)該段圓曲線參數(shù)的函數(shù)模型可以表示為

[(xi-exi)-a]2+[(yi-eyi)-b]2=r2

(1)

式中：exi和eyi為與各坐標(biāo)觀測值相對應(yīng)的隨機(jī)誤差。

NTLS的極小值條件是所有輸入非線性系統(tǒng)的觀測誤差平方和最小。采用數(shù)學(xué)語言描述可表達(dá)為

(2)

式中：ex和ey分別為與向量x和y相對應(yīng)的誤差向量。用GNSS采集道路曲線坐標(biāo)時(shí)，通常認(rèn)為觀測條件相同，屬于等精度觀測，觀測誤差服從期望為0，方差為σ2的正態(tài)分布，它們對應(yīng)的權(quán)矩陣都是單位矩陣。

2 3種極小值條件的比較

附有參數(shù)的條件平差法、最小二乘法和NTLS，都是在某一極小值條件下，估計(jì)道路圓曲線參數(shù)的最佳近似值。它們的極小值條件分別為式(3)、式(4)和式(2)：

(3)

(4)

(5)

通過觀察它們的幾何關(guān)系，可以發(fā)現(xiàn)，觀測誤差exi、eyi和半徑誤差eri，分別是直角三角形的垂直邊和斜邊，它們滿足勾股定理，即

(6)

(7)

比較式(7)和式(2)，不難發(fā)現(xiàn)它們具有一致性，這說明，同附有參數(shù)的條件平差法一樣，NTLS的極小值條件具有明確的幾何意義，比最小二乘法定義的極小值條件更合理。

3 NTLS反演道路圓曲線參數(shù)的解算

將式(1)變形，得到非線性總體最小二乘平差函數(shù)式(8)：

f(ex,ey,a,b,r)=[(x-ex)-a]2+[(y-ey)-b]2-r2

(8)

令a=a0+δa，b=b0+δb，r=r0+δr，對式(8)進(jìn)行線性化，整理后可得

AE+Bξ-W=0

(9)

式中：

A=[A1,A2]

(10)

(11)

(12)

ξ=[δaδbδr]T

(13)

(14)

用上述符號，可將NTLS平差極小值條件式(2)重新表達(dá)為

ETPE=min

(15)

式中：P=I2m，I2m是大小為2m×2m的單位矩陣。

采用拉格朗日乘數(shù)法求解NTLS，首先需要構(gòu)造目標(biāo)函數(shù)，根據(jù)式(9)和式(15)，可將該目標(biāo)函數(shù)描述為

F(E,ξ,λ)=ETPE+2λT(AE+Bξ-W)

(16)

式中：λ是拉格朗日乘子。

將式(16)分別對E，ξ和λ求偏導(dǎo)數(shù)，并令它們等于0，即

(17)

(18)

(19)

式中：帶波浪線和小帽的符號分別表示最佳預(yù)測值和最佳估計(jì)值。

對式(17)～式(19)進(jìn)行整理，可以得到參數(shù)修正值和誤差預(yù)測值的計(jì)算公式，即

(20)

(21)

驗(yàn)后單位權(quán)方差的計(jì)算公式為

(22)

4 NTLS反演圓曲線參數(shù)算法

Step1：計(jì)算參數(shù)初值。

Step2：計(jì)算參數(shù)估計(jì)值。

1)計(jì)算觀測量和半徑的估計(jì)值

3)計(jì)算參數(shù)修正值和誤差預(yù)測值

4)更新參數(shù)值

Step3：判斷是否終止迭代。

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

為驗(yàn)證、測試筆者提出的算法，設(shè)計(jì)模擬數(shù)據(jù)(實(shí)驗(yàn)1)和實(shí)測數(shù)據(jù)(實(shí)驗(yàn)2)兩個(gè)實(shí)驗(yàn)。

5.1 實(shí)驗(yàn)1

為保證模擬數(shù)據(jù)盡可能地接近實(shí)測數(shù)據(jù)，按如下規(guī)則產(chǎn)生1 000組含有隨機(jī)誤差的觀測數(shù)據(jù)。首先，設(shè)置接近實(shí)際的圓心坐標(biāo)和半徑，(a,b)=(436000,32400)，r=700 ，單位m；其次，從180°開始，在圓心角為5°所對應(yīng)的圓弧上，用參數(shù)方程，等間距計(jì)算出20個(gè)不含有誤差的圓曲線點(diǎn)坐標(biāo)；最后，給這些精確的坐標(biāo)值，附加上期望為0，中誤差為0.05 m(方差： 0.002 5 m2)的隨機(jī)誤差，模擬1 000次，得到1 000組模擬觀測數(shù)據(jù)。

分別采用重心坐標(biāo)法、最小二乘法、附有參數(shù)的條件平差法和NTLS對這1 000組數(shù)據(jù)進(jìn)行擬合，估計(jì)出圓心坐標(biāo)和半徑，將4種方法每次估計(jì)的參數(shù)與真值做差，結(jié)果如圖1。從圖1可以看出：①4種方法的估計(jì)結(jié)果存在明顯差異；②相對于其它兩種方法而言，附有參數(shù)的條件平差法和NTLS的估計(jì)結(jié)果更穩(wěn)定，特別是后者；③NTLS的估計(jì)值與真值最為接近；④附有參數(shù)的條件平差法估計(jì)的圓半徑始終存在明顯偏差。

圖1 估計(jì)值與真值之差Fig. 1 Difference between estimated value and truth value

圖2描述了附有參數(shù)的條件平差法和NTLS估計(jì)的驗(yàn)后單位權(quán)方差值，它們的平均值分別為24.379 7和0.002 7，顯然，后者更接近先驗(yàn)值。

圖2 驗(yàn)后單位權(quán)方差估計(jì)值Fig. 2 Estimated values of posteriori unit weight variance

5.2 實(shí)驗(yàn)2

實(shí)測數(shù)據(jù)實(shí)驗(yàn)，數(shù)據(jù)來源于文獻(xiàn)[3]，如表1。仍然采用上述4種方法反演道路圓曲線參數(shù)，反演結(jié)果如表2, 對應(yīng)的圖形描述如圖3。

表1 GNSS實(shí)測道路圓曲線坐標(biāo)Table 1 GNSS measured road circular curve coordinates m

表2 實(shí)驗(yàn)2參數(shù)估計(jì)結(jié)果Table 2 Parameter estimation results in experiment 2 m

圖3 實(shí)測道路圓曲線反演結(jié)果Fig. 3 Inversion results of measured road circular curve

從表2可以發(fā)現(xiàn)，4種方法的反演結(jié)果存在明顯差異，其中附有參數(shù)的條件平差法反演的圓曲線半徑明顯小于其余方法，其驗(yàn)后單位權(quán)方差遠(yuǎn)大于NTLS，這與模擬計(jì)算結(jié)果一致?？紤]到道路圓曲線的半徑一般較大，在實(shí)驗(yàn)2中采用重心坐標(biāo)法和NTLS更為合理。但是，考慮到實(shí)驗(yàn)1的結(jié)果，后者比前者不僅穩(wěn)定性更好，而且更接近真值。因此，綜合考慮，在道路圓曲線參數(shù)反演中推薦使用NTLS。

6 結(jié) 論

道路圓曲線參數(shù)的反演是恢復(fù)道路線形參數(shù)的關(guān)鍵。當(dāng)采樣點(diǎn)坐標(biāo)x和y都含有誤差時(shí)，同時(shí)考慮它們的影響，并直接對它們進(jìn)行處理，有助于提高參數(shù)估計(jì)精度，增加平差理論的嚴(yán)密性。根據(jù)前文的討論，在道路圓曲線參數(shù)反演中，可以得出如下結(jié)論：

1)NTLS的極小值條件等價(jià)于點(diǎn)到圓曲線距離平方和最小條件。

2)NTLS反演道路圓曲線參數(shù)，與其余3種方法相比較，不僅穩(wěn)定性好，而且精度更高。

3)綜合考慮，就筆者討論的4種方法而言，推薦使用NTLS反演道路圓曲線參數(shù)。

事實(shí)上，不僅觀測數(shù)據(jù)本身含有粗差，而且一些被錯(cuò)誤識別的點(diǎn)也屬于粗差，因此采用穩(wěn)健總體最小二乘估計(jì)圓曲線參數(shù)其理論更嚴(yán)密，這將是我們進(jìn)一步研究的方向。