葉玲潔，顏遠(yuǎn)青

(1.廣州市城市規(guī)劃勘測(cè)設(shè)計(jì)研究院，廣東 廣州 510060； 2.珠海大橫琴科技發(fā)展有限公司，廣東 珠海 519000)

1 引 言

在建筑物的三維重建中，屋頂面提取是其中最關(guān)鍵的部分[1]，后續(xù)工作都在屋頂面準(zhǔn)確提取的基礎(chǔ)上進(jìn)行。按照屋頂面的平整度，可以將屋頂面分為平面式屋頂和曲面式屋頂，現(xiàn)代建筑中尤其是城市建筑物，一般都較為規(guī)則，屋頂面多為平面式，而對(duì)于曲面式屋頂，難以直接用參數(shù)定義其曲面形狀，常規(guī)的做法是對(duì)其構(gòu)建Delaunay三角網(wǎng)，以三角面片的形式來擬合曲面。本文的研究重點(diǎn)主要是針對(duì)平面式屋頂面，因此對(duì)曲面式的屋頂面不進(jìn)行研究討論。

在屋頂面分割算法中，目前常用的主要是RANSAC(Random Sample Consensus，隨機(jī)抽樣一致性)算法[2]、區(qū)域增長(zhǎng)算法[3]、三維霍夫變換算法[4]、聚類算法等。區(qū)域增長(zhǎng)算法首先選定一個(gè)種子區(qū)域，根據(jù)屋頂面點(diǎn)云與種子區(qū)域的法向量夾角、空間位置關(guān)系不斷地拓展面片，該方法分割效果良好，但是算法中參數(shù)選取較為困難，同時(shí)由于法向量的計(jì)算易受噪聲點(diǎn)影響，進(jìn)而使得算法增長(zhǎng)出錯(cuò)誤的平面。RANSAC是一種隨機(jī)參數(shù)估計(jì)算法，首先從點(diǎn)云中隨機(jī)選取三個(gè)點(diǎn)，計(jì)算出平面參數(shù)，然后計(jì)算其他點(diǎn)與該平面的偏差，當(dāng)偏差小于閾值，則該點(diǎn)為局內(nèi)點(diǎn)，不斷迭代找到局內(nèi)點(diǎn)最多的平面，則為最優(yōu)平面模型。RANSAC能魯棒的估計(jì)模型參數(shù)，但是其迭代次數(shù)必須足夠多才能得到準(zhǔn)確的結(jié)果，效率低下。

在進(jìn)行平面擬合時(shí)，常用的方法是最小二乘法進(jìn)行平面擬合，但最小二乘法對(duì)于自變量與因變量無法區(qū)分的情況，適應(yīng)性差[5]，在實(shí)際中，一般認(rèn)為接近垂直于地面的屋頂面是不存在的，以Z坐標(biāo)作為因變量來對(duì)屋頂面進(jìn)行最小二乘擬合平面。同時(shí)，在屋頂面點(diǎn)云存在噪聲點(diǎn)的情況下，對(duì)擬合結(jié)果影響很大。

本文利用PCA(Principal Component Analysis，主成分分析法)來計(jì)算點(diǎn)云的法向量，并以區(qū)域增長(zhǎng)法來對(duì)平面進(jìn)行分割。PCA算法可以有效減小噪聲點(diǎn)對(duì)法向量計(jì)算的影響，改善區(qū)域增長(zhǎng)的結(jié)果。

2 平面擬合

本文中采用PCA的方法對(duì)點(diǎn)云進(jìn)行平面擬合。PCA是一種數(shù)學(xué)變換的方法，利用降維的思想在變換中保持變量的總方差不變，將給定的一組變量線性變換為另一組不相關(guān)的變量，并且使變換后的第一變量的方差最大，即第一主成分，其他分量的方差依次遞減[6]。在本文中的變量為三維點(diǎn)坐標(biāo)的集合，其變量為X、Y、Z三個(gè)坐標(biāo)值，則經(jīng)過變換后，應(yīng)有三個(gè)主成分，對(duì)于一個(gè)空間平面，在平行于平面的方向上點(diǎn)集分布最為離散，方差最大，在垂直于平面的方向上，點(diǎn)集分布最為集中，方差最小，即空間平面的第三主成分為垂直于空間平面的向量。由于平面擬合最關(guān)鍵的為法向量的擬合，利用PCA得到點(diǎn)集的第三主成分，即能進(jìn)一步擬合出平面方程，如圖1所示。

圖1 PCA變換原理

對(duì)于在坐標(biāo)系XYZ下的待擬合平面點(diǎn)云，利用主成分分析法對(duì)其進(jìn)行分析，可得到三個(gè)按照從大到小排列的特征值λ1、λ2、λ3，對(duì)應(yīng)的主分量分別為V1、V2、V3，其中V1和V2組成了待擬合平面的一組基，V3與V1、V2正交，為垂直于待擬合平面的法向量。如圖1，在XYZ坐標(biāo)系下的點(diǎn)云，經(jīng)過主成分分析后，三個(gè)主成分分量V1、V2、V3組成了新坐標(biāo)系X′Y′Z′的三個(gè)基，V1和V2為平面X′O′Z′的一組基，V3為O′Z′方向的基，即所擬合平面的法向量。

PCA過程如下：

(1)特征中心化。即每一維的數(shù)據(jù)都減去該維的均值，變換之后每一維的均值都變成了零。特征中心化后的點(diǎn)集P，如式(1)，其中，xi、yi、zi為中心化后點(diǎn)坐標(biāo)：

(1)

(2)計(jì)算三個(gè)坐標(biāo)的協(xié)方差矩陣。協(xié)方差矩陣C為：

(2)

其中，cov(x，y)為x坐標(biāo)和y坐標(biāo)的協(xié)方差，cov(x，x)為x坐標(biāo)的方差，協(xié)方差計(jì)算公式如式(3)，xi、yi為中心化后點(diǎn)坐標(biāo)：

(3)

當(dāng)協(xié)方差大于零時(shí)說明x和y是正相關(guān)關(guān)系，協(xié)方差小于零時(shí)x和y是負(fù)相關(guān)關(guān)系，協(xié)方差為零時(shí)x和y相互獨(dú)立。

(3)計(jì)算協(xié)方差矩陣C的特征值和特征向量。所計(jì)算出來的特征值按照從大到小排序，分別為λ1、λ2、λ3，其所對(duì)應(yīng)的特征向量分別為ξ1、ξ2、ξ3。顯然，兩個(gè)較大λ所對(duì)應(yīng)的特征向量ξ1、ξ2為待擬合平面的一組基，而ξ3為待擬合平面的法向量，其三個(gè)分量分別為a、b、c。

若已知待擬合平面經(jīng)過點(diǎn)p(x0，y0，z0)p(x0，y0，z0)，則擬合平面為式(4)：

a(x-x0)+b(y-y0)+c(z-z0)=0

(4)

采用主成分分析法擬合平面的方法，對(duì)于存在噪聲點(diǎn)的情況，也能很好的擬合出結(jié)果。因?yàn)樵谝粋€(gè)平面點(diǎn)云中，噪聲點(diǎn)偏離平面的距離相對(duì)于平面的范圍較小，對(duì)擬合結(jié)果的影響可以忽略。

3 法向量計(jì)算

傳統(tǒng)的平面分割算法中，對(duì)于點(diǎn)的法向量計(jì)算通常是采用以點(diǎn)為圓心，r為半徑的一個(gè)鄰域，取鄰域內(nèi)的所有點(diǎn)云進(jìn)行最小二乘擬合，這種方法在平滑變化的曲面上，計(jì)算的結(jié)果不會(huì)有太大偏差，但是在邊緣處，法向量難以確定，其結(jié)果是不準(zhǔn)確的[7]。

本文中，首先對(duì)點(diǎn)的法向量進(jìn)行計(jì)算。在選擇點(diǎn)云的鄰域點(diǎn)時(shí)，把當(dāng)前表面情況考慮進(jìn)去，即對(duì)于在同一個(gè)平面上，且空間上相鄰近，則認(rèn)為是鄰域點(diǎn)，若不在同一個(gè)平面上，即使空間上相鄰近，也不視為鄰域點(diǎn)[8]。

圖2 初始鄰域點(diǎn)的選擇

法向量計(jì)算的主要步驟如下：

(1)定義所求點(diǎn)的r鄰域，搜索到其所有鄰域點(diǎn)，并進(jìn)行平面擬合，擬合結(jié)果為Φinit。其中，如圖2，初始的鄰域點(diǎn)為在以所求點(diǎn)為球心，r為半徑的區(qū)域內(nèi)的所有點(diǎn)，都認(rèn)為是所求點(diǎn)的鄰域點(diǎn)；

(2)計(jì)算鄰域點(diǎn)到Φinit的距離di的中誤差σ，將距離di小于兩倍中誤差的點(diǎn)權(quán)值設(shè)為1，距離di大于兩倍中誤差的點(diǎn)權(quán)值設(shè)為0。重新擬合平面Φinit；權(quán)值計(jì)算公式如(5)：

(5)

(3)多次擬合，直到所擬合的平面參數(shù)不再改變或者改變量小于閾值，此時(shí)擬合的平面為Φi。則Φi的法向量為所求點(diǎn)的法向量。

圖3 法向量求解過程

法向量的求解過程，也就是所求點(diǎn)處平面的擬合過程，可以用圖3表示，初始擬合的平面采用了r鄰域內(nèi)的所有點(diǎn)云進(jìn)行擬合，其結(jié)果為Φinit，鄰域內(nèi)與Φinit的距離小于2σ的權(quán)值為1，如圖中黃色點(diǎn)，將權(quán)值為1的點(diǎn)重新擬合，多次擬合，最終結(jié)果為Φi，從圖3中可看出，此時(shí)所采用的用來擬合Φi的鄰域點(diǎn)為待求點(diǎn)的r鄰域內(nèi)，且在同一個(gè)平面內(nèi)的點(diǎn)。如圖4，將以待求點(diǎn)為球心，r為半徑的區(qū)域內(nèi)，與待求點(diǎn)在同一個(gè)平面內(nèi)的點(diǎn)，作為待求點(diǎn)的鄰域點(diǎn)，以最終確定的鄰域點(diǎn)來計(jì)算所求點(diǎn)的法向量。

圖4 最后鄰域點(diǎn)的確定

4 平面分割

本文中平面分割的主要算法為區(qū)域增長(zhǎng)，主要步驟如下[9]：

(1)生成分割結(jié)果列表L，種子點(diǎn)隊(duì)列S，當(dāng)前區(qū)域點(diǎn)列表Q；

(2)計(jì)算所有點(diǎn)的法向量和曲率。點(diǎn)pi對(duì)應(yīng)的法向量、曲率分別為ni、ci；

(3)選擇曲率最小且未聚類點(diǎn)作為種子點(diǎn)，將其添加到種子點(diǎn)隊(duì)列S和當(dāng)前平面點(diǎn)列表Q中；

④若種子點(diǎn)隊(duì)列S為空時(shí)，當(dāng)前平面提取完成，將當(dāng)前平面點(diǎn)列表Q保存到分割結(jié)果列表L中；重置Q。

(5)最后的分割結(jié)果保存在列表L中。

在將點(diǎn)云數(shù)據(jù)進(jìn)行平面分割后，同一個(gè)平面上的點(diǎn)云將被劃分到同一個(gè)聚類中。此時(shí)，對(duì)其進(jìn)行平面擬合，計(jì)算出平面參數(shù)。做法是采用PCA算法計(jì)算得到平面的法向量，即平面參數(shù)a、b、c，將點(diǎn)集內(nèi)所有點(diǎn)求其中心，以其中心作為平面上點(diǎn)，計(jì)算出平面參數(shù)d。

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

圖5 人字形屋頂面的區(qū)域增長(zhǎng)算法分割結(jié)果

圖6 錐形屋頂面的區(qū)域增長(zhǎng)算法分割結(jié)果

圖7 L形屋頂面的區(qū)域增長(zhǎng)算法分割結(jié)果

圖5、圖6、圖7，為三組不同數(shù)據(jù)，采用傳統(tǒng)的區(qū)域增長(zhǎng)對(duì)數(shù)據(jù)進(jìn)行平面分割的結(jié)果，其中，angle為角度閾值，可以看出，在平面上無噪聲或者噪聲較少的情況下，區(qū)域增長(zhǎng)算法提取出來的平面內(nèi)部點(diǎn)較為連續(xù)，不會(huì)出現(xiàn)空洞，但在邊緣處的點(diǎn)云無法很好地分割開。這是因?yàn)檫吘壧幍狞c(diǎn)云法向量計(jì)算困難，無法得到其正確的法向量，因此邊緣處的點(diǎn)云分割較為零碎，而本文中設(shè)定，當(dāng)點(diǎn)數(shù)少于閾值時(shí)，不判別為一個(gè)平面，即一個(gè)平面上，必須有足夠多的點(diǎn)，才能被識(shí)別出。本文設(shè)置閾值為20個(gè)點(diǎn)。同時(shí)，區(qū)域增長(zhǎng)的參數(shù)不容易設(shè)置，對(duì)于不同形狀的點(diǎn)云需要不同的參數(shù)，使得建筑物平面點(diǎn)云自動(dòng)化分割難以實(shí)現(xiàn)，制約了建筑物自動(dòng)化三維重建的可能。

圖8 錐形屋頂面的RANSAC算法分割結(jié)果

圖9 L形屋頂面的RANSAC算法分割結(jié)果

圖8、圖9，為兩組數(shù)據(jù)，采用了RANSAC算法對(duì)數(shù)據(jù)進(jìn)行平面分割的結(jié)果，其中，T為閾值?？梢钥闯觯琑ANSAC在參數(shù)合適的時(shí)候能夠較好地分割出平面。但是，RANSAC提取出的每個(gè)平面并不是等價(jià)的，第一個(gè)被提取出的平面總是會(huì)擁有更多的點(diǎn)，這一點(diǎn)在圖8的T=1.0 m時(shí)表現(xiàn)比較明顯。此外，參數(shù)不合適時(shí)，出現(xiàn)各種錯(cuò)誤情況，如T=0.1 m時(shí)，平面中出現(xiàn)很多空洞，這是由于機(jī)載LiDAR掃描的點(diǎn)云在高程方向有一定的波動(dòng)，但被RANSAC錯(cuò)誤地濾掉；在第一組T=1.5 m和第二組T=1.0 m時(shí)，則出現(xiàn)全局最優(yōu)，RANSAC只把全局內(nèi)，分配到當(dāng)前平面點(diǎn)數(shù)最多的判別為平面。

圖10本文算法結(jié)果與影像對(duì)比

圖10為本文算法分割結(jié)果與實(shí)際影像對(duì)比，可以看出分割的結(jié)果較為準(zhǔn)確，在平面相交的邊緣處也沒有因?yàn)辄c(diǎn)云法向的不確定性而導(dǎo)致的錯(cuò)誤分割或者邊緣缺失的情況。

而對(duì)于圖11類型的穹頂建筑，本文算法無法適應(yīng)。

圖11 穹頂屋頂建筑

6 結(jié) 語

本文簡(jiǎn)單介紹了常用的平面分割算法，包括三維霍夫變換、RANSAC算法和區(qū)域增長(zhǎng)算法。三維霍夫變換需要將點(diǎn)云變換到屬性空間，計(jì)算量較大，一般較少使用；RANSAC算法需多次迭代，效率低下，在參數(shù)不合適的時(shí)候，也容易出現(xiàn)全局最優(yōu)而導(dǎo)致分割錯(cuò)誤，此外，若第一個(gè)提取的平面擁有更多的點(diǎn)也導(dǎo)致后續(xù)提取的平面質(zhì)量不斷下降；區(qū)域增長(zhǎng)算法對(duì)邊緣處的法向量計(jì)算比較困難，隨之出現(xiàn)的是邊緣處的點(diǎn)云無法很好地分割到平面上。

本文采用PCA算法擬合局部平面的方法，并根據(jù)點(diǎn)到平面的距離設(shè)置權(quán)值，迭代多次后得到最準(zhǔn)確的平面及其所對(duì)應(yīng)的法向量，由于所選擇的鄰域一般較小，所迭代次數(shù)相對(duì)于RANSAC算法大大減少。但在屋頂面點(diǎn)云數(shù)量較少的情況下，噪聲點(diǎn)對(duì)屋頂面法向量計(jì)算影響變大，容易導(dǎo)致法向量計(jì)算出錯(cuò)，從而使得分割結(jié)果錯(cuò)誤。計(jì)算出法向量后，再以區(qū)域增長(zhǎng)算法來對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行分割，經(jīng)過實(shí)驗(yàn)，其結(jié)果準(zhǔn)確，分割效果良好。根據(jù)本文算法進(jìn)行點(diǎn)云平面分割，其結(jié)果可應(yīng)用于建筑物的建模、點(diǎn)云分類等。本文算法主要針對(duì)平面點(diǎn)云的分割，對(duì)于平頂、坡頂甚至多平面復(fù)合型建筑物屋頂點(diǎn)云分割效果較好，但對(duì)于穹頂、曲面建筑物屋頂點(diǎn)云無法適用，筆者將繼續(xù)深入研究并加以改進(jìn)。

[1] 胡偉，盧小平，李珵等. 基于改進(jìn)RANSAC算法的屋頂激光點(diǎn)云面片分割方法[J]. 測(cè)繪通報(bào)，2012(11)：31～34.

[2] 李娜，馬一薇，楊洋等. 利用RANSAC算法對(duì)建筑物立面進(jìn)行點(diǎn)云分割[J]. 測(cè)繪科學(xué)，2011，36(5)：144～145.

[3] 李寶順，岑紅燕，包亞萍等. 基于平面提取的點(diǎn)云數(shù)據(jù)分割算法[J]. 計(jì)算機(jī)應(yīng)用與軟件，2014，31(7)：145～148.

[4] 黃先鋒. 利用機(jī)載LIDAR數(shù)據(jù)重建3D建筑物模型的關(guān)鍵技術(shù)研究[D]. 武漢：武漢大學(xué)，2006.

[5] 浮丹丹，周紹光，徐洋等. 基于主成分分析的點(diǎn)云平面擬合技術(shù)研究[J]. 測(cè)繪工程，2014，23(4)：20～23.

[6] 王松，夏紹瑋. 一種魯棒主成分分析(PCA)算法[J]. 系統(tǒng)工程理論與實(shí)踐，1998(1)：10～14.

[7] 楊斌. 機(jī)載LiDAR點(diǎn)元數(shù)據(jù)建筑物半自動(dòng)提取方法研究[D]. 阜新：遼寧工程技術(shù)大學(xué)，2012.

[8] Kim C，Habib A，Pyeon M，et al. Segmentation of Planar Surfaces from Laser Scanning Data Using the Magnitude of Normal Position Vector for Adaptive Neighborhoods[J]. Sensors，2016，16(2)：140.

[9] 楊洋，張永生，張皓等. 基于LiDAR數(shù)據(jù)的建筑物自動(dòng)化重建[J]. 測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2010，27(2)：123～126.

[10] 姜如波. 基于三維激光掃描技術(shù)的建筑物模型重建[J]. 城市勘測(cè)，2013，3：113～116.