蔣榮華，劉 超

(安徽理工大學(xué)測(cè)繪學(xué)院，安徽 淮南 232000)

一種三維激光掃描點(diǎn)云擬合的抗差加權(quán)整體最小二乘法

蔣榮華，劉 超

(安徽理工大學(xué)測(cè)繪學(xué)院，安徽 淮南 232000)

針對(duì)三維激光掃描中點(diǎn)云不等精度且易受粗差影響的問(wèn)題，提出了一種基于入射角定權(quán)的抗差加權(quán)總體最小二乘的擬合方法。該方法在采用入射角定權(quán)的基礎(chǔ)上，進(jìn)行基于標(biāo)準(zhǔn)化殘差和中位數(shù)的抗差加權(quán)整體最小二乘估計(jì)，獲得待定參數(shù)估值，并通過(guò)Gauss-Newton迭代算法，推導(dǎo)了模型的迭代計(jì)算方法。以平面擬合和球面擬合為例，分別通過(guò)仿真數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)對(duì)算法進(jìn)行驗(yàn)證，結(jié)果表明，對(duì)于含有粗差的點(diǎn)云，新方法可以獲得更為理想的參數(shù)估值，其性能優(yōu)于抗差整體最小二乘和加權(quán)整體最小二乘，可以更好地進(jìn)行三維激光掃描的點(diǎn)云擬合。

三維激光掃描；點(diǎn)云；加權(quán)整體最小二乘；抗差估計(jì)；入射角；定權(quán)

Abstract： In point clouds fitting of terrestrial laser scanner，the data observed with unequal and the accuracy is affected greatly by gross error．In order to overcome this shortcoming，a robust weighted total least squares method based on the theory of incidence angle weighting is adopted in this paper，and an incidence angle cosine formula to spherical point cloud is deduced．So as to get the coefficient and inspection，this method uses the theory of incident angle to determine weights and utilizes the standardized residuals to construct the weight factor function，and the square root of the variance component estimator with robustness is obtained by the median method．There use experiments based on the simulation and measured data of plane fitting and the sphere fitting as example to verification algorithm，the experiment indicates that the robust WTLS(RWTLS) method exhibits satisfactory robustness，the accuracy of the obtained parameters is high，it is superior to the robust total least squares and weighted total least squares method．It has perfect performance on point clouds of terrestrial laser scanner．

Keywords： terrestrial laser scanning；point clouds；weighted total least squares；robust estimation；incidence angle;weighing

三維激光掃描作為一種較新的數(shù)據(jù)獲取技術(shù)，具有高精度、高效率、無(wú)接觸、快速等優(yōu)勢(shì)，已成為數(shù)據(jù)獲取的重要手段。其中如何利用點(diǎn)云數(shù)據(jù)建立精確三維模型是其應(yīng)用的一個(gè)關(guān)鍵問(wèn)題。國(guó)內(nèi)外學(xué)者對(duì)點(diǎn)云建模進(jìn)行了大量研究，并取得一些有益成果。點(diǎn)云擬合最初采用最小二乘法(LS)，考慮系數(shù)矩陣含有誤差，Colub H等提出了整體最小二乘方法[1](total least squares，TLS)，后Schaffirin發(fā)表了帶有等式約束的TLS方法[2]；但整體最小二乘方法沒(méi)有考慮數(shù)據(jù)不等精度問(wèn)題，為解決此問(wèn)題，Schaffrin提出在整體最小二乘基礎(chǔ)上加入權(quán)陣的加權(quán)整體最小二乘迭代方法(weighted total least squares，WTLS)[3]，較好地解決了變量不等精度的整體最小二乘問(wèn)題。在三維激光掃描領(lǐng)域，后續(xù)研究者也進(jìn)行了深入研究[4-7]。陳瑋嫻提出了根據(jù)點(diǎn)云激光反射強(qiáng)度及對(duì)系數(shù)陣A列向量部分修正引入權(quán)陣的加權(quán)整體最小二乘法，很好地解決了平面靶標(biāo)和球面靶標(biāo)擬合問(wèn)題[8]；蒼桂華提出了入射角定權(quán)的加權(quán)總體最小二乘平面擬合方法[9]，獲得較高擬合精度。但這些方法均沒(méi)有考慮觀測(cè)數(shù)據(jù)含有異常值情況，因此，需采用更為穩(wěn)健的參數(shù)估計(jì)方法。龔循強(qiáng)提出了基于IGG權(quán)函數(shù)的穩(wěn)健加權(quán)總體最小二乘方法，獲得了更可靠、更優(yōu)參數(shù)估值[10]；王彬提出了抗差加權(quán)整體最小二乘的牛頓-高斯算法，此算法利用標(biāo)準(zhǔn)化殘差進(jìn)行檢驗(yàn)量的構(gòu)建，并利用中位數(shù)法獲取單位權(quán)中誤差估值，具有很好抗差性[11]。然而，以上抗差方法目前鮮有應(yīng)用到三維激光領(lǐng)域，因此，針對(duì)三維激光掃描點(diǎn)云數(shù)據(jù)擬合中數(shù)據(jù)不等精度及含有異常值的問(wèn)題，在結(jié)合抗差整體最小二乘及入射角定權(quán)方法基礎(chǔ)上，本文提出基于三維激光掃描點(diǎn)云數(shù)據(jù)的抗差加權(quán)整體最小二乘方法。此方法采用點(diǎn)云數(shù)據(jù)入射角余弦值cosθ作為擬合的權(quán)重，利用標(biāo)準(zhǔn)化殘差構(gòu)造權(quán)因子函數(shù)，利用中位數(shù)法獲得單位權(quán)中誤差估值，并應(yīng)用于三維激光掃描平面、耙球的擬合，經(jīng)實(shí)例證明此方法可以有效解決點(diǎn)云數(shù)據(jù)含有粗差及不等精度擬合問(wèn)題。

1 抗差加權(quán)整體最小二乘

1.1 函數(shù)模型及基本原理

由于地面三維激光現(xiàn)實(shí)掃描場(chǎng)景中含有大量平面特征及激光靶標(biāo)多為平面球面，因此本文以平面和球面為例。

對(duì)于空間平面點(diǎn)云數(shù)據(jù)，方程式為

z=ax+by+c

(1)

式中，a、b、c為待求的平面擬合參數(shù)；x、y、z為平面點(diǎn)云數(shù)據(jù)三維坐標(biāo)。

將平面方程寫成矩陣形式

Z=AX

(2)

式中

(3)

對(duì)于球面，方程為

(4)

式中，a、b、c為球心坐標(biāo)；x、y、z為球面點(diǎn)云數(shù)據(jù)三維坐標(biāo)；r為球半徑。球面方程展開(kāi)式為

x2+y2+z2=2xa+2yb+2zc+r2-a2-b2-c2

(5)

將x2+y2+z2、2x、2y、2z看作新的觀測(cè)值，a、b、c、r2-a2-b2-c2看作新的參數(shù)，可寫成矩陣形式

Z=AX

(6)

式中

(7)

同時(shí)顧及觀測(cè)向量Z和系數(shù)矩陣A的誤差，建立EIV模型

Z-eZ=A-EAX

(8)

相應(yīng)的隨機(jī)模型為

(9)

(10)

式中

(11)

其中權(quán)因子函數(shù)Rii計(jì)算公式如下

(12)

(13)

(14)

根據(jù)協(xié)因數(shù)傳播率得

(15)

1.2 權(quán)陣的確定

在三維激光掃描點(diǎn)云擬合中，由于點(diǎn)云數(shù)據(jù)是不等精度觀測(cè)，本文根據(jù)各點(diǎn)精度確定權(quán)值，相關(guān)文獻(xiàn)提出點(diǎn)云數(shù)據(jù)入射角越小，其點(diǎn)位精度越高[6-7]，其權(quán)值越大，因此本文引入入射角余弦值作為權(quán)重。平面點(diǎn)云數(shù)據(jù)觀測(cè)點(diǎn)權(quán)值為

(16)

式中，Pi為觀測(cè)點(diǎn)權(quán)值；θi為i點(diǎn)入射角。

(17)

式中，(a，b，-1)為擬合平面法向量；xi、yi、zi含義同式(4)。同理，推出球面點(diǎn)云觀測(cè)點(diǎn)入射角余弦值為

(18)

式中，(xi-a，yi-b，zi-c)為擬合球面法向量；xi、yi、zi含義同式(4)。

1.3Gauss-Neuton迭代算法

設(shè)第i次迭代計(jì)算后，所得參數(shù)估值為X(i)，預(yù)測(cè)殘差矩陣為EA(i)。

由于EIV模型是非線性的，將式(8)在(X(i),EA(i))處線性化展開(kāi)，保留一階導(dǎo)數(shù)項(xiàng)，過(guò)程如下

Z-eZ=(A-EA)X= (A-EA(i))X(i)+(A-EA(i))dX-dEAX(i)=AX(i)+(A-EA(i))dX-(EA(i)+dEA)X(i)=AX(i)+A(i)dX-EAX(i)

(19)

轉(zhuǎn)換并移項(xiàng)得

(20)

結(jié)合EIV模型估計(jì)準(zhǔn)則，構(gòu)造拉格朗日目標(biāo)函數(shù)

(21)

根據(jù)Eular-Lagrange必要條件對(duì)變量eZ、eA、dX、K求導(dǎo)并令其等于零，有

(22)

(23)

(24)

(25)

代入式(22)得

(26)

(27)

迭代步驟如下：

(2) 根據(jù)式(17)、式(18)分別計(jì)算平面和球面各點(diǎn)余弦值，組成相關(guān)權(quán)陣及協(xié)因數(shù)陣。

(3) 從i=0開(kāi)始依次計(jì)算

A(i)=A-EA(i)

(28)

(29)

(30)

(31)

(32)

(33)

(5) 將WTLS迭代獲得的參數(shù)估值和殘差向量作為初值，將第k-1次迭代獲得的參數(shù)估值和殘差向量，以及先驗(yàn)協(xié)因數(shù)陣QZ、QA代入式(15)更新相應(yīng)殘差協(xié)因數(shù)陣，再將殘差向量和殘差協(xié)因數(shù)陣代入式(13)和式(14)獲得第k次迭代標(biāo)準(zhǔn)化殘差。

(6) 利用式(12)和式(13)更新等價(jià)協(xié)因數(shù)陣，從而獲得第k+1次參數(shù)估值和殘差向量。

2 算例分析

2.1 仿真數(shù)據(jù)試驗(yàn)

為驗(yàn)證算法可行性，首先用仿真數(shù)據(jù)進(jìn)行試驗(yàn)。以平面為例，設(shè)擬合平面方程為z=x+2y+1，給定參數(shù)真值為：a=1，b=2，c=1。在Matlab中在坐標(biāo)(xi,yi,zi)上加入均值為u=0，方差為σ2=1/100×cosθi的隨機(jī)誤差(其中i表示點(diǎn)號(hào)，θi表示入射角余弦值)，從中任取選取30個(gè)點(diǎn)，得到一組模擬點(diǎn)云數(shù)據(jù)。

在模擬數(shù)據(jù)中加入位置隨機(jī)產(chǎn)生的粗差，其絕對(duì)值一般介于標(biāo)準(zhǔn)差的10～30倍之間。設(shè)計(jì)方案如下：加入粗差前分別進(jìn)行加權(quán)整體最小二乘(WTLS)和抗差加權(quán)整體最小二乘(RWTLS)估計(jì)，然后分別加入1、2、3個(gè)粗差，比較兩種方法參數(shù)估計(jì)結(jié)果，見(jiàn)表1。

表1 平面模擬數(shù)據(jù)結(jié)果

表1中max(da)、max(db)、max(dc)分別為da、db、dc的最大擬合參數(shù)估值與真值差。

由表1可看出：當(dāng)沒(méi)有粗差時(shí)，RWTLS方法與WTLS方法精度相同，均獲得較高精度參數(shù)解；當(dāng)含粗差時(shí)，RWTLS方法擬合參數(shù)估值與真值差明顯小于WTLS，且與不含粗差時(shí)精度相當(dāng)，參數(shù)估值準(zhǔn)確度較高。

為直觀比較兩種方法的參數(shù)擬合準(zhǔn)確性，在粗差為3時(shí)進(jìn)行100次模擬，比較100次模擬兩種方法da、db差值序列，結(jié)果如圖1所示。

圖1 粗差個(gè)數(shù)為3時(shí)WTLS和RWTLS參數(shù)差值序列

由圖1可看出，RWTLS獲得的參數(shù)差值較WTLS小、更接近真值。

2.2 實(shí)測(cè)數(shù)據(jù)試驗(yàn)

為進(jìn)一步驗(yàn)證算法準(zhǔn)確性，采用實(shí)測(cè)數(shù)據(jù)并加入抗差整體最小二乘方法(RTLS)進(jìn)行計(jì)算。試驗(yàn)采用中海達(dá)LS300三維激光掃描儀對(duì)某校測(cè)繪學(xué)院紅樓和布設(shè)球靶標(biāo)(半徑為0.072 5 m)的學(xué)校羅馬柱廣場(chǎng)進(jìn)行掃描，得到平面和球面點(diǎn)云數(shù)據(jù)，從中選取平面和球面數(shù)據(jù)，手動(dòng)刪除大量冗余點(diǎn)，得到如圖2點(diǎn)云數(shù)據(jù)，各選取50個(gè)點(diǎn)進(jìn)行試驗(yàn)，方案同上，參數(shù)擬合結(jié)果見(jiàn)表2、表3。

圖2 平面和球面靶標(biāo)點(diǎn)云數(shù)據(jù)

粗差個(gè)數(shù)算法^a^b^c^σ20(×10-5)0WTLS72.8899-40.0484629.20310.1337RTLS72.5053-39.8360625.91470.3437RWTLS72.8899-40.0484629.20310.13371WTLS74.8093-41.1344645.38140.5142RTLS72.8102-40.0222628.33650.2317RWTLS72.8167-40.0057628.59910.13142WTLS76.3429-41.9158659.18282.2635RTLS72.8038-40.0254628.27630.2106RWTLS72.8635-40.0572628.73330.15473WTLS82.6095-45.2382714.45192.9177RTLS73.3441-40.3188632.90490.3640RWTLS73.0012-40.1253630.00640.1675

表3 球面實(shí)測(cè)數(shù)據(jù)計(jì)算結(jié)果

表中dr2為半徑估值平方與真值差值。由表2、表3可看出，對(duì)于平面和球面擬合，當(dāng)不含粗差時(shí)，RWTLS與WTLS方法精度高于RTLS方法，說(shuō)明采用入射角定權(quán)方法可獲得較高精度參數(shù)解；當(dāng)含粗差時(shí)，WTLS擬合精度明顯下降，參數(shù)解波動(dòng)較大，且隨著粗差個(gè)數(shù)的增加擬合精度越來(lái)越低，而RTLS方法和RWTLS方法擬合精度與不含粗差時(shí)相當(dāng)，單位權(quán)中誤差變化較小，優(yōu)于WTLS方法，且RWTLS方法獲得單位權(quán)中誤差和半徑平方差值比RTLS方法小，精度優(yōu)于RTLS方法。

3 結(jié) 語(yǔ)

針對(duì)三維激光掃描數(shù)據(jù)常規(guī)擬合方法精度受粗差影響較大情況，本文結(jié)合抗差整體最小二乘和入射角定權(quán)方法提出一種以入射角定權(quán)的抗差加權(quán)整體最小二乘方法，解決了三維激光掃描點(diǎn)云數(shù)據(jù)擬合問(wèn)題。理論推導(dǎo)和試驗(yàn)分析結(jié)果表明：

(1) 當(dāng)數(shù)據(jù)不含粗差時(shí)，本文采用的RWTLS方法擬合精度高于RTLS方法，可以得出RWTLS方法采用的入射角定權(quán)方法可以提高擬合參數(shù)精度，從而獲得較優(yōu)參數(shù)解。

(2) 當(dāng)數(shù)據(jù)含粗差時(shí)，WTLS方法擬合精度明顯下降，且隨著粗差個(gè)數(shù)的增加，WTLS擬合精度越來(lái)越低，而RWTLS和RTLS擬合精度幾乎不變且與不含粗差時(shí)精度相當(dāng)。

綜上，以入射角定權(quán)的RWTLS具有較好抗差性和較高精度參數(shù)解，能夠解決三維激光掃描點(diǎn)云擬合中數(shù)據(jù)含有異常值及不等精度問(wèn)題。

[1] GOLUB H，AN LOAN CH．An Analysis of the Total Least Squares Problem[J]．SIAM Journal on Numerical Analysis，1980，17(6)：883-893．[2] SCHAFFRIN B，US Y A．On Total Least-square Adjustment with Constraints[J]．IAG-Symp，2005，128：417-421．

[3] SCHAFFRIN B，IESER A．On Weighted Total Least-square Adjustment for Liner Regression[J]．Journal of Geodesy，2008，82(7)：415-421．

[4] 王樂(lè)洋，許才軍，魯鐵定．病態(tài)加權(quán)總體最小二乘平差的嶺估計(jì)解法[J]．武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2011，35(11)：1346-1350．

[5] AMIRI S，JAZAERI S．Weighted Total Least Squares Formulated by Standard Least Squares Theory[J]．Journal of Geodesy，2012，2(2)：113-124．

[6] VIAHID M．On Weighted Total Least-squares for Geodetic Transformations[J]．Journal of Geodesy，2012，86 (5)：359-367．

[7] JAZAERI S，AMIRI S，SHARIFI M A．Iterative Algorithm for Weighted Total Least Squares Adjustment[J]．Survey Review，2014，46(334)：19-27．

[8] 陳瑋嫻，陳義，袁慶，等．加權(quán)總體最小二乘在三維激光標(biāo)靶擬合中的應(yīng)用[J]．大地測(cè)量與地球動(dòng)力學(xué)，2010，30(5)：90-96．

[9] 蒼桂華，李明峰，岳建平．以入射角定權(quán)的點(diǎn)云數(shù)據(jù)加權(quán)總體最小二乘平面擬合研究[J]．大地測(cè)量與地球動(dòng)力學(xué)，2014，34(3)：95-98．

[10] 龔循強(qiáng)，李志林．穩(wěn)健加權(quán)總體最小二乘法[J]．測(cè)繪學(xué)報(bào)，2014，43(9)：888-894．

[11] 王彬，李建成，高井祥，等．抗差加權(quán)整體最小二乘模型的牛頓-高斯算法[J]．測(cè)繪學(xué)報(bào)，2015，44(6)：602-607．

[12] 陳瑋嫻，袁慶．抗差總體最小二乘[J]．大地測(cè)量與地球動(dòng)力學(xué)，2012，32(6)：111-114．

[13] 官云蘭，劉紹堂，周世健，等．基于整體最小二乘的穩(wěn)健點(diǎn)云數(shù)據(jù)平面擬合[J]．大地測(cè)量與地球動(dòng)力學(xué)，2011，31(5)：80-83．

[14] 楊元喜.抗差估計(jì)的概念和任務(wù)[J].測(cè)繪通報(bào)，2004(4)：36-39.

[15] 楊元喜，宋力杰，徐天河．大地測(cè)量相關(guān)觀測(cè)抗差估計(jì)理論[J]．測(cè)繪學(xué)報(bào)，2002，31(2)：95-99．

[16] 袁慶，樓立志，陳瑋嫻．基于加權(quán)整體最小二乘的點(diǎn)云數(shù)據(jù)平面擬合法[J]．測(cè)繪通報(bào)，2011(3)：1-3．

RobustWeightedTotalLeastSquaresMethodforTerrestrialLaserScanningPointCloudFitting

JIANG Ronghua，LIU Chao

(Anhui University of Science and Technology，Huainan 232000，China)

P237

A

0494-0911(2017)09-0037-05

蔣榮華，劉超.一種三維激光掃描點(diǎn)云擬合的抗差加權(quán)整體最小二乘法[J].測(cè)繪通報(bào)，2017(9)：37-41.

10.13474/j.cnki.11-2246.2017.0283.

2017-05-18

國(guó)家自然科學(xué)基金(41404004)；安徽省博士后基金(2015B044)；安徽理工大學(xué)科研啟動(dòng)基金(11152)

蔣榮華(1994—)，女，碩士生，研究方向?yàn)槿S激光掃描點(diǎn)云擬合及抗差估計(jì)。E-mail：844510389@qq．com