基于加權總體最小二乘法的點云平面擬合

蒼桂華，岳建平

引 言

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

1 IWTLS方法及其函數模型

1.1 EIV 模型

設點云數據所建立3維空間平面方程式形式為:

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

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

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

根據系數矩陣特點引入權陣P0，PX和PY。P0是3×3矩陣，代表系數矩陣A的列向量權陣;PX是n×n矩陣，代表系數矩陣A的行向量權陣;PY是n×n矩陣，代表向量 Y的權陣。P0，PX，PY相對應的協因素矩陣為 Q0，QX，QY，即:

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

式中，?表示“kronecker積”。

IWTLS估算準則為:

1.2 IWTLS權陣的設置

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

設點云在x，y和z3個方向等精度獲取，對于平面的系數陣列向量和觀測向量中，有σx=σy=σz。結合系數矩陣A的特點，設置相應權陣。

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

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

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

(1)根據(8)式計算出各點強度值，并根據(9)式設置相關矩陣 P0，PX，PY。

2 算例及分析

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

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

將單位權中誤差、點到擬合面的最大距離以及平面擬合精度作為估算方法優劣的評判指標。標準反射板、普通木板、水泥模板3種樣本數據采用不同平面擬合方法得到的相關結果如表1所示。

Fig.1 Plane fitting of point clouds in experiments

Table 1 Results of plane fitting samples

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

3 小結

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

［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).