許輝熙,薛萬蓉2,程　熙3

(1. 四川建筑職業(yè)技術(shù)學(xué)院測(cè)量工程研究所, 四川 德陽 618000; 2. 四川建筑職業(yè)技術(shù)學(xué)院

測(cè)繪工程系, 四川 德陽 618000; 3. 成都理工大學(xué),四川 成都 610059)

XUHuixi,XUEWanrong,CHENGXi

一種擬合三維空間直線的新方法

許輝熙1,2,3,薛萬蓉2,程熙3

(1. 四川建筑職業(yè)技術(shù)學(xué)院測(cè)量工程研究所, 四川 德陽 618000; 2. 四川建筑職業(yè)技術(shù)學(xué)院

測(cè)繪工程系, 四川 德陽 618000; 3. 成都理工大學(xué),四川 成都 610059)

ANewFittingMethodofThree-dimensionalSpatialStraightLines

XUHuixi,XUEWanrong,CHENGXi

摘要：提出了一種基于穩(wěn)健總體最小二乘的三維空間直線擬合新方法。該方法以加權(quán)總體最小二乘為基礎(chǔ)，通過刪除點(diǎn)到擬合直線距離過大的觀測(cè)點(diǎn)來抵抗粗差的影響，獲得空間直線參數(shù)的穩(wěn)健估計(jì)值。試驗(yàn)表明，當(dāng)觀測(cè)數(shù)據(jù)中不含有誤差時(shí)，加權(quán)最小二乘法、加權(quán)總體最小二乘法和本文方法的參數(shù)估計(jì)值高度一致；當(dāng)觀測(cè)數(shù)據(jù)包含粗差時(shí)，本方法的參數(shù)估計(jì)值明顯更接近真實(shí)值。

引文格式： 許輝熙,薛萬蓉,程熙. 一種擬合三維空間直線的新方法[J].測(cè)繪通報(bào)，2015(9)：28-31.DOI:10.13474/j.cnki.11-2246.2015.0271

關(guān)鍵詞：穩(wěn)健估計(jì);加權(quán)總體最小二乘;三維空間直線擬合

中圖分類號(hào)：P258

文獻(xiàn)標(biāo)識(shí)碼：B

文章編號(hào)：0494-0911(2015)09-0028-04

收稿日期：2015-03-30

基金項(xiàng)目：國家自然科學(xué)基金(41301488)；四川省教育廳自然科學(xué)項(xiàng)目(15ZB0448;15ZB0452)

作者簡介：許輝熙(1979—),男,博士,副教授,主要從事地理學(xué)、空間信息技術(shù)應(yīng)用研究。E-mail：529382949@qq.com

一、引言

三維空間直線擬合是工程測(cè)量和工業(yè)測(cè)量等工程應(yīng)用領(lǐng)域的常見問題。三維空間直線可以用點(diǎn)向式表達(dá)為

(1)

式中，(x0,y0,z0)為空間直線通過的已知點(diǎn)；(f,g,h)為空間直線的方向矢量。擬合平面直線是典型的線性問題，而擬合空間直線是一個(gè)明顯的非線性問題，不能簡單直接地借用擬合平面直線的加權(quán)最小二乘法(leastsquares，LS)和加權(quán)總體最小二乘法[1](weightedtotalleastsquares，WTLS)。

要運(yùn)用LS法進(jìn)行空間直線擬合，需要先將空間直線描繪成兩個(gè)空間平面的交線形式,或?qū)⑵浯怪蓖队暗阶鴺?biāo)平面上，然后擬合平面或平面直線，最后重建空間直線[2-3]。也可以先擬合通過該空間直線的任何一個(gè)平面，然后在該平面上擬合平面直線，最后將平面直線還原為空間直線[4]。很顯然，這些方法含有某一個(gè)坐標(biāo)分量不含誤差的假設(shè)，不符合實(shí)際情況。文獻(xiàn)[1]為消除上述假設(shè)，采用總體最小二乘法(totalleastsquares，TLS)進(jìn)行空間直線擬合，取得了良好效果，但未針對(duì)含有粗差的情況進(jìn)行討論。三維空間直線擬合問題是一個(gè)非線性問題，因此可以采用非線性最小二乘法(nonlinearleastsquares,NLS)。文獻(xiàn)[5]提出采用高斯-牛頓迭代算法擬合空間直線，這些迭代算法對(duì)初值有一定的要求，并且可能出現(xiàn)發(fā)散，從而得不到正確的參數(shù)估計(jì)結(jié)果。文獻(xiàn)[6]討論了穩(wěn)健總體最小二乘法(robusttotalleastsquares,RTLS)在三維相似坐標(biāo)變換中的應(yīng)用。

目前，關(guān)于用RTLS法進(jìn)行空間直線擬合的討論并不多見。本文提出采用非線性拉格朗日函數(shù)法推導(dǎo)基于WTLS的一種擇優(yōu)錄取RTLS算法。算法首先以含有粗差的觀測(cè)數(shù)據(jù)擬合直線，然后計(jì)算各觀測(cè)點(diǎn)到直線的距離，擇優(yōu)錄取距離較小的點(diǎn)重新擬合空間直線。

二、總體最小二乘擬合三維空間直線

將式(1)進(jìn)行簡單變形，表達(dá)成如下的兩個(gè)空間平面[1]

(2)

式中，ξ1=f/h;ξ2=x0-(fz0)/h；ξ3=g/h；ξ4=y0-(gz0)/h。因此，只要擬合出平面的4個(gè)參數(shù)就可以重建六參數(shù)形式的空間直線對(duì)稱方程。假設(shè)有t(t>2)個(gè)實(shí)測(cè)空間點(diǎn)數(shù)據(jù)，用這些觀測(cè)點(diǎn)擬合平面參數(shù)的矩陣形式為

(3)

根據(jù)式(3)可以得到WTLS擬合空間直線的函數(shù)模型，即

L-eL=(B-EB)ξ

(4)

式中，eL和EB分別為觀測(cè)值的隨機(jī)誤差矢量和矩陣，大小分別為m×1和m×4(m=2t)。從式(4)中可以明顯發(fā)現(xiàn)，此時(shí)系數(shù)矩陣具有很強(qiáng)的結(jié)構(gòu)特性。如果采用SVD分解解算將會(huì)使得系數(shù)矩陣中的常數(shù)變量分配到不應(yīng)有的誤差改正值，不同位置的同一變量得不到相同的改正值，這與理論不符，因此解算理論不嚴(yán)密。為了解決這種結(jié)構(gòu)性問題，將式(4)進(jìn)行改寫，使誤差矩陣只與z坐標(biāo)分量誤差有關(guān)，即

L-eL=Bξ-(ξT?Im)Hez

(5)

式中，Hez=eB=vec(EB)。H為常數(shù)矩陣，vec()表示矩陣?yán)庇?jì)算，從左往右，將后一列加到前一列的后面。方程中的誤差滿足如下的隨機(jī)屬性

(6)

(7)

式中，Qx=Qy=Qz=Q，Q表示觀測(cè)點(diǎn)之間的協(xié)因數(shù)矩陣。根據(jù)式(3)—式(7)的描述，可以得到如下的非線性拉格朗日函數(shù)

Φ=eLPleL+ezPzez+

2λT(L-eL-Bξ+(ξT?Im)Hez)

(8)

要使上述非線性拉格朗日函數(shù)的極小值存在，則Φ關(guān)于eL、ez、λ和ξ的偏導(dǎo)數(shù)等于零，即

(9)

對(duì)式(9)進(jìn)行推算整理，可以得到如下的參數(shù)估計(jì)公式

(10)

式中

并且可以得到誤差改正項(xiàng)的計(jì)算公式，即

(11)

(12)

根據(jù)式(3)和式(11)可以得到x和y坐標(biāo)分量的誤差改正值，即

(13)

很顯然，誤差改正量等于觀測(cè)點(diǎn)到擬合直線各坐標(biāo)分量上的距離值。因此，點(diǎn)到直線的距離即是各坐標(biāo)分量改正值的平方和的根，即

(14)

如此則可以根據(jù)點(diǎn)的觀測(cè)精度確定極限距離，將式(14)確定的距離與極限距離相比，找出在極限距離范圍內(nèi)的點(diǎn)。最后的加權(quán)單位權(quán)方差估計(jì)值計(jì)算公式為

(15)

式中，q表示擇優(yōu)錄取點(diǎn)的總數(shù)量。WTLS問題是一個(gè)簡單的非線性問題，單位權(quán)方差估計(jì)值一般是有偏的。

三、穩(wěn)健總體最小二乘擬合算法

隨著三維激光掃描儀、激光跟蹤儀及其他先進(jìn)儀器在工程測(cè)量中的廣泛使用，空間直線點(diǎn)云數(shù)據(jù)的獲取成為可能。在點(diǎn)的采集過程中由于受到外界環(huán)境的影響，點(diǎn)云中可能包含有異常點(diǎn)。在利用WTLS法進(jìn)行空間直線的參數(shù)擬合時(shí)沒有考慮這些異常點(diǎn)的影響，擬合的空間直線參數(shù)與目標(biāo)空間直線參數(shù)間存在較大差異。而RTLS法可以消除這些異常點(diǎn)，使得擬合的空間直線更接近目標(biāo)空間直線。其具體計(jì)算流程如下：

1) 根據(jù)給定的坐標(biāo)觀測(cè)值和精度確定平差需要的觀測(cè)矢量、系數(shù)矩陣及協(xié)方差因子。

2) 假設(shè)EB=0，采用LS法計(jì)算參數(shù)的初始值。

3) 采用式(10)循環(huán)計(jì)算參數(shù)的估計(jì)值直至相鄰兩次估計(jì)值之差的2范數(shù)小于指定限差，本文的限差設(shè)置為1Ε-10。

4) 采用式(11)—式(12)計(jì)算殘差矢量，并將殘差矢量描述成各坐標(biāo)分量的殘差矢量。

5) 采用式(14)計(jì)算觀測(cè)點(diǎn)到擬合直線的距離。

6) 采用式(15)計(jì)算加權(quán)單位權(quán)方差估計(jì)值。

8) 利用擇優(yōu)錄取的點(diǎn)重復(fù)步驟1)—步驟7)；直至相鄰兩次錄取的點(diǎn)相同。

9) 輸出估計(jì)的4個(gè)參數(shù)，將參數(shù)代入式(2)中恢復(fù)出空間直線的對(duì)稱方程。

四、試驗(yàn)分析

在一條已知的空間直線

(16)

上隨機(jī)產(chǎn)生200個(gè)空間點(diǎn)，如圖1所示。將該空間

圖1　原始模擬數(shù)據(jù)

直線投影到xoz和yoz坐標(biāo)平面上，空間直線可以通過這兩個(gè)垂直投影面表達(dá)，即

(17)

為了驗(yàn)證本文方法的可行性，設(shè)計(jì)如下兩個(gè)試驗(yàn)方案。

1) 方案1：在不含有人為粗差的原始模擬數(shù)據(jù)中，對(duì)x、y、z坐標(biāo)分量附加上期望為零、中誤差為0.03的隨機(jī)誤差，產(chǎn)生一組模擬觀測(cè)值；分別采用LS、WTLS和RTLS法估計(jì)式(2)的參數(shù)。

2) 方案2：在原始模擬數(shù)據(jù)中隨機(jī)抽出10個(gè)點(diǎn)，在這些點(diǎn)上加上人為粗差，同樣給這些粗差點(diǎn)和其余原始模擬點(diǎn)附加上期望為零、中誤差為0.03的隨機(jī)誤差，并采用上述3種方法估計(jì)模型參數(shù)。

對(duì)方案1和方案2計(jì)算1000次，表1和表2分別列出了方案1和方案2的參數(shù)估計(jì)平均值和方差估計(jì)平均值。圖2分別表達(dá)了方案1和方案2的參數(shù)估計(jì)值與真值的差值。方案1和方案2每次的方差估計(jì)值分別描繪如圖3所示。

表1　LS、WTLS和RTLS估計(jì)結(jié)果與真值的比較(方案1)

表2　LS、WTLS和RTLS估計(jì)結(jié)果與真值的比較(方案2)

從表1和圖2(a)中不難發(fā)現(xiàn)，在不含有粗差的情況下，3種方的參數(shù)估計(jì)值非常接近，并且與真值的差異幾乎相同。但是圖3(a)則說明RTLS法得到的方差估計(jì)值比LS和WTLS的小，在不含有人為粗差的情況下，一些隨機(jī)誤差較大的點(diǎn)被RTLS法剔除了。

表2和圖2(b)說明在含有粗差的情況下RTLS法估計(jì)的參數(shù)更接近真實(shí)值。從表2的最后一列和圖3(b)可以發(fā)現(xiàn)，此時(shí)3種方法得到的方差估計(jì)值存在較大的差異，同方案1一樣，RTLS計(jì)算的方差估計(jì)值最小，這說明RTLS法能夠有效的剔除粗差。

五、結(jié)束語

在將空間直線擬合問題轉(zhuǎn)換成擬合兩個(gè)相交平面問題后，LS估計(jì)沒有顧及z坐標(biāo)軸的觀測(cè)誤差，估計(jì)理論不夠嚴(yán)密。雖然WTLS法能夠解決上述問題，也能夠處理系數(shù)矩陣中的結(jié)構(gòu)性問題，但是對(duì)于含有粗差的情況卻無能為力。為了解決這些問題，本文提出了采用RTLS法估計(jì)空間直線參數(shù)，用非線性拉格朗日乘數(shù)法推導(dǎo)了RTLS的計(jì)算公式，設(shè)計(jì)了一種可行的RTLS算法。模擬試驗(yàn)結(jié)果證明了本文算法能夠有效剔除粗差，提高了參數(shù)和驗(yàn)后方差的估計(jì)精度。

圖2　LS、WTLS和RTLS估計(jì)參數(shù)與真值之差

圖3　LS、WTLS和RTLS的方差估計(jì)值

參考文獻(xiàn)：

[1]姚宜斌,黃書華,孔建,等.空間直線擬合的整體最小二乘算法[J]. 武漢大學(xué)學(xué)報(bào)：信息科學(xué)版,2014,39(5): 571-574.

[2]襲楊.空間直線擬合的一種方法[J].齊齊哈爾大學(xué)學(xué)報(bào): 自然科學(xué)版,2009, 25(2): 64-68.

[3]王偉鋒, 溫耐.空間直線擬合研究[J].許昌學(xué)院學(xué)報(bào), 2010,29(5): 37-39.

[4]郭際明,向巍,尹洪斌.空間直線擬合的無迭代算法[J]. 測(cè)繪通報(bào), 2011(2): 24-26.

[5]王解先,季凱敏.工業(yè)測(cè)量擬合[M].北京:測(cè)繪出版社,2008.

[6]陳義,陸玨.以三維坐標(biāo)轉(zhuǎn)換為例解算穩(wěn)健總體最小二乘方法[J]. 測(cè)繪學(xué)報(bào), 2012, 41(5): 715-722.

[7]李博峰,沈云中.基于等效殘差積探測(cè)粗差的方差-協(xié)方差分量估計(jì)[J]. 測(cè)繪學(xué)報(bào)，2011, 40(1): 10-14.

[8]欒學(xué)晨，楊必勝，張?jiān)品?基于結(jié)構(gòu)模式的道路網(wǎng)結(jié)點(diǎn)匹配方法[J]. 測(cè)繪學(xué)報(bào)， 2013, 42(4): 608-614.