基于加權(quán)總體最小二乘法的點(diǎn)云平面擬合

蒼桂華，岳建平

引 言

點(diǎn)云數(shù)據(jù)平面擬合是地面3維激光掃描技術(shù)應(yīng)用中常見問題，如建筑物點(diǎn)云數(shù)據(jù)的特征線提取等［1］。最小二乘(least squares，LS)法是常用的平面擬合方法。設(shè)利用點(diǎn)云數(shù)據(jù)建立的平面方程的矩陣形式為Y=AX，最小二乘法是假設(shè)誤差e只存在于觀測向量Y中，建立經(jīng)典的高斯-馬爾科夫(Gaussian-Markov，G-M)模型對誤差方程式進(jìn)行求解，獲得平面參量估值。然而由于模型誤差、人為誤差、儀器誤差等因素影響點(diǎn)云數(shù)據(jù)在觀測數(shù)據(jù)x，y和z這3個方向上均存在誤差，使得包含觀測數(shù)據(jù)的系數(shù)矩陣A也含有誤差。因此，利用最小二乘法進(jìn)行平面數(shù)據(jù)擬合的結(jié)果并不是最優(yōu)，而是有偏的［2］。針對這種觀測向量和系數(shù)矩陣均包含誤差的模型，即所謂的變量中的誤差(error-in-variables，EIV)模型，GOLUB等人提出了總體最小二乘(total least squares，TLS)估計(jì)方法［3］，但總體最小二乘估計(jì)僅在設(shè)計(jì)矩陣和殘差元素均服從獨(dú)立等精度分布時才是最優(yōu)估計(jì)［4］。實(shí)際點(diǎn)云數(shù)據(jù)中各點(diǎn)坐標(biāo)精度是不等的，因此簡單的總體最小二乘方法并非最優(yōu)估計(jì)［5］。為了解決矩陣和殘差的不等精度估計(jì)問題，MARKOVSKY等人提出了加權(quán)總體最小二乘(weighted totalleastsquares，WTLS)方 法［6］。SCHAFFRIN等人則進(jìn)一步擴(kuò)展了WTLS方法，詳細(xì)介紹了相關(guān)權(quán)陣的設(shè)計(jì)方法以及算法的步驟［7］。本文中在SCHAFFRIN等人提出的WTLS方法基礎(chǔ)上，根據(jù)點(diǎn)云數(shù)據(jù)中各點(diǎn)反射強(qiáng)度值與其點(diǎn)位精度關(guān)系，確定各點(diǎn)的平面擬合權(quán)值，得到強(qiáng)度加權(quán)總體最小二乘(intensity weighted total least squares，IWTLS)的平面擬合方法，并通過均質(zhì)性不同的3種樣本點(diǎn)云數(shù)據(jù)，對該方法的適用性進(jìn)行研究。

1 IWTLS方法及其函數(shù)模型

1.1 EIV 模型

設(shè)點(diǎn)云數(shù)據(jù)所建立3維空間平面方程式形式為:

式中，a，b和c為待求的平面擬合參量。

將(1)式寫成矩陣形式為:

如果同時考慮觀測向量Y和系數(shù)矩陣A中的誤差，則建立EIV函數(shù)模型［7］:

根據(jù)系數(shù)矩陣特點(diǎn)引入權(quán)陣P0，PX和PY。P0是3×3矩陣，代表系數(shù)矩陣A的列向量權(quán)陣;PX是n×n矩陣，代表系數(shù)矩陣A的行向量權(quán)陣;PY是n×n矩陣，代表向量 Y的權(quán)陣。P0，PX，PY相對應(yīng)的協(xié)因素矩陣為 Q0，QX，QY，即:

由 Q0，QX可以得到 QA，PA:

式中，?表示“kronecker積”。

IWTLS估算準(zhǔn)則為:

1.2 IWTLS權(quán)陣的設(shè)置

激光反射強(qiáng)度值與入射角關(guān)系為Ii=I0cosαi(i=1，2，…，n)，Ii代表入射角為 αi時點(diǎn)的強(qiáng)度值，I0為垂直入射時(αi=0)點(diǎn)的反射強(qiáng)度值。入射角越小，點(diǎn)位精度越高，點(diǎn)的反射強(qiáng)度值越大［8-10］。因此強(qiáng)度值越大，參與擬合的權(quán)重應(yīng)越大。本次實(shí)驗(yàn)數(shù)據(jù)為.PTS格式，以12bit記錄強(qiáng)度值，其強(qiáng)度值范圍為［－2047，2048］。設(shè)記錄的原始強(qiáng)度值為Ii'，按照下式將其值變?yōu)椋?，1］之間，構(gòu)成各點(diǎn)的擬合權(quán)值:

設(shè)點(diǎn)云在x，y和z3個方向等精度獲取，對于平面的系數(shù)陣列向量和觀測向量中，有σx=σy=σz。結(jié)合系數(shù)矩陣A的特點(diǎn)，設(shè)置相應(yīng)權(quán)陣。

式中，P0的第3個對角元素為0，表示系數(shù)矩陣A中的第3列不需要改正，其余對角線元素為1，表示系數(shù)矩陣A中的第1列和第2列的數(shù)據(jù)列中的元素是等精度獲取的;PX，PY與強(qiáng)度值有關(guān)。

1.3 平面參量求解及精度評定［7，9-10］

目前解決WTLS問題主要采用基于拉格朗日乘數(shù)法的迭代解算方法［7］，計(jì)算步驟如下。

(1)根據(jù)(8)式計(jì)算出各點(diǎn)強(qiáng)度值，并根據(jù)(9)式設(shè)置相關(guān)矩陣 P0，PX，PY。

2 算例及分析

為了驗(yàn)證IWTLS方法的適用性，采用了均質(zhì)性不同的3種平面樣本進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)中的平面樣本分別為標(biāo)準(zhǔn)反射板(反射率為90%)、普通木板及一般的水泥建材模板。3種平面樣本中，標(biāo)準(zhǔn)反射板均質(zhì)性最好;普通木板次之;建材水泥模板均質(zhì)性最差。利用徠卡ScanStation2 C10分別對樣本進(jìn)行掃描，獲取點(diǎn)云數(shù)據(jù)(如圖1所示)。

根據(jù)點(diǎn)云數(shù)據(jù)特點(diǎn)確定平面方程式形式［11］。分別利用LS法、TLS法和IWTLS法對各個樣本點(diǎn)云數(shù)據(jù)進(jìn)行平面擬合，獲得平面擬合參量，，以及單位權(quán)中誤差。設(shè)擬合平面上點(diǎn)的個數(shù)為n，計(jì)算出各點(diǎn)i(i=1，2，…，n)到擬合面的距離di，獲得點(diǎn)到擬合面的最大距離di，max，根據(jù)下式計(jì)算出平面擬合精度:

將單位權(quán)中誤差、點(diǎn)到擬合面的最大距離以及平面擬合精度作為估算方法優(yōu)劣的評判指標(biāo)。標(biāo)準(zhǔn)反射板、普通木板、水泥模板3種樣本數(shù)據(jù)采用不同平面擬合方法得到的相關(guān)結(jié)果如表1所示。

Fig.1 Plane fitting of point clouds in experiments

Table 1 Results of plane fitting samples

從表1可以看出:對于均質(zhì)性較好的標(biāo)準(zhǔn)反射率板和普通木板，利用IWTLS方法獲得的3個精度判定指標(biāo)(σ^0，σ^p，di，max)值均要比 LS 方法和 TLS 方法的相應(yīng)結(jié)果小得多。以普通木頭為例，利用IWTLS方法得到的單位權(quán)中誤差比LS方法和TLS方法分別提高了99%和25%，平面擬合精度比LS方法和TLS方法分別提高了63%和40%，點(diǎn)到擬合面的最大距離也由LS方法和TLS方法的8.3mm和5.8mm降為5.1mm。然而對于一般的建材水泥模板，由于均質(zhì)性較差，IWTLS 方法計(jì)算出的σ^0，σ^p和di，max3 個精度評判指標(biāo)分別為 0.000323m，3.2mm和6.5mm，雖好于LS方法的相應(yīng)結(jié)果0.037647m，4.2mm和 8.5mm，卻比 TLS方法的相關(guān)結(jié)果0.000304m，1.8mm 和5.8mm 差。這是由于此時各點(diǎn)強(qiáng)度值的差異更多由于材質(zhì)不同造成，其強(qiáng)度值已經(jīng)不能代表其點(diǎn)位精度，因此，利用IWTLS方法獲得的相關(guān)結(jié)果差于總體最小二乘方法。

3 小結(jié)

點(diǎn)云數(shù)據(jù)平面擬合中WTLS法雖然從理論上較LS法和TLS法合理，但在擬合時應(yīng)注意各點(diǎn)擬合權(quán)值的設(shè)置。確定的各點(diǎn)擬合權(quán)值應(yīng)與其點(diǎn)位精度一致，即點(diǎn)位精度高，擬合權(quán)值應(yīng)越大。如果擬合權(quán)重與其實(shí)際點(diǎn)位精度情況不一致，會直接影響WTLS的效果。本文中的強(qiáng)度加權(quán)總體最小二乘法對于均質(zhì)性較好的點(diǎn)云平面效果明顯，擬合精度較高，而對于均質(zhì)性較差的點(diǎn)云平面效果不佳，此時應(yīng)采用TLS方法進(jìn)行平面擬合。

［1］ YU H X，WU K，AO J F，et al.Extraction of building’s feature lines based on 3-D laser scanning technology［J］.Laser Technology，2012，36(4):553-556(in Chinese).

［2］ QIU W N，TAO B Z，YAO Y B，et al.The theory and method of surveying data processing［M］.Wuhan:Wuhan University Press，2008:161-175(in Chinese).

［3］ GOLUB G H，van LOAN C F.An analysis of the total least squares problem［J］.SIAM Journal on Numerical Analysis，1980，17(6):883-893.

［4］ van HUFFEL S，VANDEWALLE J.The total least squares problem:computational aspects and analysis［M］.Philadelphia，USA:Society for Industrial and Applied Mathematics，1991:263-283.

［5］ ZHOU Y J，DENG C H.Weighted and unweighted total least square methods and applications to heteroscedastic 3-D coordinate transformation［J］.Geomatics and Information Science of Wuhan University，2012，37(8):976-979(in Chinese).

［6］ MARKOVSKY I，RASTELLO M，PREMOLI A，et al.The element-wise weighted total least-squares problem［J］.Computational Statistics and Data Analysis，2006，50(1):181-209.

［7］ SCHAFFRIN B，WIESER A.On weighted total least-square adjustment for linear regression ［J］.Journal of Geodesy，2008，82(7):415-421.

［8］ ZHANG Y.Research on point cloud processing of terrestrial laser scanning［D］.Wuhan:Wuhan University，2008:43-47(in Chinese).

［9］ CHEN W X，CHEN Y，YUAN Q，et al.Application of weighted total least squares to target fitting of three-dimensional laser scanning［J］.Journal of Geodesy and Geodynamics，2010，30(5):90-96(in Chinese).

［10］ YUAN Q，LOU L Z，CHEN W X.Appling weight total leastsquares to the plane point cloud fitting of terrestrial laser scanning［J］.Bulletin of Surveying and Mapping，2011(3):1-3(in Chinese).

［11］ WANG J X，JI K M.Industrial surveying fitting［M］.Beijing:Surveying and Mapping Press，2007:43-45(in Chinese).