鄒紀(jì)偉, 劉德兒, 康 翔, 蘭小機(jī), 楊 鵬, 劉靖鈺

(江西理工大學(xué) 土木與測(cè)繪工程學(xué)院, 江西 贛州 341000)

0 引 言

三維激光掃描技術(shù)自問世以來已被應(yīng)用于各個(gè)領(lǐng)域, 如文物修復(fù)、建筑物的三維建模、智慧城市的建設(shè)和大面積的地形測(cè)量等[1]。地面激光掃描儀能夠快速高精度地獲取地物的表面信息, 相比于機(jī)載激光掃描儀和航空攝影測(cè)量, 其獲取的數(shù)據(jù)更加精細(xì)[2], 由于點(diǎn)云數(shù)據(jù)量龐大, 若直接對(duì)其進(jìn)行操作, 工作量較大且耗時(shí)長, 因此需對(duì)源數(shù)據(jù)進(jìn)一步處理, 點(diǎn)云分割作為數(shù)據(jù)處理的關(guān)鍵環(huán)節(jié), 其分割精度直接影響到后續(xù)點(diǎn)云三維重建、特征提取、地物識(shí)別等操作。

常見的點(diǎn)云分割算法有RANSAC算法、歐氏聚類算法、區(qū)域增長算法、最小割算法等。RANSAC算法作為隨機(jī)參數(shù)的估計(jì)算法[3], 由Fishier等于1981年提出, 該算法在直線擬合、平面擬合、計(jì)算圖形方面應(yīng)用廣泛, 在點(diǎn)云分割中, 該算法通常用于提取平面點(diǎn)云：李娜等[4]將RANSAC算法運(yùn)用于點(diǎn)云的分割處理, 提取出了建筑物的特征立面；賀亦峰等[5]則通過設(shè)置不同閾值, 使用RANSAC算法將建筑立面中的墻面、窗戶、陽臺(tái)等要素單獨(dú)分割開來, 再通過歐氏聚類完成獨(dú)立個(gè)體的分割?；跉W氏聚類的點(diǎn)云分割是根據(jù)不同點(diǎn)間的歐氏距離來判斷是否屬于同類地物點(diǎn)：吳燕雄等[6]使用歐氏聚類加上平滑閾值的方法對(duì)地物實(shí)現(xiàn)分類提取, 能夠高效快速地達(dá)到分割效果, 但若地物種類過多, 閾值則較難設(shè)置, 易出現(xiàn)過分割或欠分割現(xiàn)象。區(qū)域增長算法[7-10]則是根據(jù)法向量及曲率實(shí)現(xiàn)不同地物分割：王競(jìng)雪等[11]提出了一種對(duì)機(jī)載LiDAR數(shù)據(jù)中的建筑物提取, 通過區(qū)域生長算法與主成分分析相結(jié)合的方法實(shí)現(xiàn)高精度提取；李仁忠等[12]提出了通過將估算的曲率最小點(diǎn)設(shè)置為種子點(diǎn)的改進(jìn)區(qū)域生長算法, 使分割的精確性及可靠性得到提高；孫金虎等[13]提出了基于最小生成樹(MST)改進(jìn)的區(qū)域生長算法, 能夠得到平滑的分割邊界。

上述常見的點(diǎn)云分割算法, 若單獨(dú)應(yīng)用于大型室外場(chǎng)景點(diǎn)云, 都存在一定的局限性, RANSAC算法適用于建筑立面細(xì)節(jié)特征的分割提取, 倘若建筑物結(jié)構(gòu)較為復(fù)雜, 則需設(shè)置不同的平面模型參數(shù)進(jìn)行分割, 操作較為繁瑣。歐氏聚類算法適用于相距較遠(yuǎn)的不同地物點(diǎn)云分割, 若不同地物點(diǎn)云出現(xiàn)粘連問題, 則很難將其單獨(dú)分割開, 閾值設(shè)置恰當(dāng)與否, 直接影響點(diǎn)云分割效果。區(qū)域增長算法因其根據(jù)法向量及曲率對(duì)散亂點(diǎn)云進(jìn)行分類, 因此對(duì)樹、桿狀物分割效果良好。針對(duì)激光點(diǎn)云中的典型地物, 找到合適的融合分割算法對(duì)其進(jìn)行分割, 將有效解決過分割或欠分割問題。

綜上, 本文提出基于區(qū)域增長與歐氏聚類相結(jié)合的點(diǎn)云分割算法：首先對(duì)數(shù)據(jù)預(yù)處理, 完成噪點(diǎn)去除及數(shù)據(jù)精簡, 采用漸近式形態(tài)學(xué)濾波去除地面點(diǎn), 再通過區(qū)域增長算法將樹、桿狀物和建筑類地物單獨(dú)分割開, 最后用歐氏聚類實(shí)現(xiàn)對(duì)典型地物的有效分割, 解決傳統(tǒng)點(diǎn)云分割算法易出現(xiàn)過分割或欠分割問題, 為后續(xù)點(diǎn)云數(shù)據(jù)處理提供基礎(chǔ)。

1 數(shù)據(jù)預(yù)處理

通過地面激光掃描儀采集到的數(shù)據(jù)由于受環(huán)境因素的影響, 不可避免地產(chǎn)生噪聲, 影響后續(xù)實(shí)驗(yàn)的分割效果, 因此本文在分割前先對(duì)源數(shù)據(jù)進(jìn)行去噪處理。具體地，統(tǒng)計(jì)濾波通過對(duì)每個(gè)點(diǎn)的鄰域進(jìn)行統(tǒng)計(jì)分析, 計(jì)算該點(diǎn)至其鄰域點(diǎn)的平均距離, 得到高斯分布圖, 將平均距離到標(biāo)準(zhǔn)范圍之外的點(diǎn)定義為離群點(diǎn), 使用時(shí)僅需設(shè)置K近鄰搜索個(gè)數(shù)和標(biāo)準(zhǔn)差倍數(shù)就可將噪聲點(diǎn)從數(shù)據(jù)中去除。本文采集數(shù)據(jù)的噪聲符合這一特性, 因此采用該算法對(duì)噪聲點(diǎn)進(jìn)行有效移除。經(jīng)濾波后, 由于數(shù)據(jù)仍然過密, 影響后續(xù)分割算法的運(yùn)行效率, 在保證數(shù)據(jù)特征信息的前提下, 采用八叉樹對(duì)數(shù)據(jù)進(jìn)行精簡。該算法首先對(duì)點(diǎn)云創(chuàng)建一個(gè)八叉樹, 在樹中搜索并記錄真實(shí)的葉節(jié)點(diǎn)數(shù)據(jù), 然后將節(jié)點(diǎn)數(shù)和對(duì)應(yīng)的數(shù)據(jù)作為平均聚類算法K值和初始聚類中心, 最終計(jì)算出每個(gè)點(diǎn)的均方根曲率、所有點(diǎn)的均方根曲率均值、每個(gè)點(diǎn)到聚類中心的歐氏距離以及每個(gè)點(diǎn)到聚類中心的歐氏距離的平均值, 根據(jù)計(jì)算結(jié)果精簡點(diǎn)云數(shù)據(jù)。

本文實(shí)驗(yàn)數(shù)據(jù)通過RIEGL VZ-1000掃描儀獲取, 其中源數(shù)據(jù)點(diǎn)個(gè)數(shù)為1 133 991, 通過對(duì)數(shù)據(jù)進(jìn)行去噪及精簡, 同時(shí)手動(dòng)對(duì)無效點(diǎn)進(jìn)行移除, 最終得到的結(jié)果見圖1, 其中統(tǒng)計(jì)濾波K近鄰搜索個(gè)數(shù)設(shè)置為6, 標(biāo)準(zhǔn)差倍數(shù)設(shè)置為1.0, 而精簡點(diǎn)云細(xì)分級(jí)別設(shè)置成10, 去噪后點(diǎn)個(gè)數(shù)為1 115 968(圖1b)，精簡及手動(dòng)去除無效點(diǎn)后點(diǎn)個(gè)數(shù)為109 130(圖1c), 可見預(yù)處理后的點(diǎn)云噪點(diǎn)被較好地移除, 地物更加清晰, 便于分割處理, 且數(shù)據(jù)精簡后將更好地提高算法的運(yùn)行效率。

圖1 數(shù)據(jù)預(yù)處理

2 點(diǎn)云分割算法

2.1 地面點(diǎn)移除

由于采集數(shù)據(jù)時(shí)不可避免夾雜地面點(diǎn)云, 為了確保點(diǎn)云分割的質(zhì)量, 需要對(duì)該區(qū)域點(diǎn)云進(jìn)行有效去除。本文研究數(shù)據(jù)為室外大型場(chǎng)景點(diǎn)云, 雖然地形較為平坦, 但仍夾雜地形起伏較大區(qū)域, 因而采用漸近式形態(tài)學(xué)濾波去除地面點(diǎn)。

形態(tài)學(xué)濾波的主要操作為腐蝕與膨脹, 該濾波算法在二維圖像中運(yùn)用較為廣泛, 其中先腐蝕后膨脹稱為開運(yùn)算, 該操作可對(duì)小區(qū)域噪聲進(jìn)行去除, 同時(shí)可對(duì)圖像進(jìn)行平滑處理; 而先膨脹后腐蝕稱為閉運(yùn)算, 該操作對(duì)圖像細(xì)小空洞有較強(qiáng)的修復(fù)能力。在激光點(diǎn)云中, 通過腐蝕運(yùn)算可獲得指定區(qū)域點(diǎn)的最小高程值, 而膨脹操作則獲得該區(qū)域內(nèi)點(diǎn)的最大高程值。借助不斷增長的濾波窗口, 通過開運(yùn)算操作, 漸近式形態(tài)學(xué)濾波即可實(shí)現(xiàn)地面點(diǎn)與非地面點(diǎn)的分離, 該濾波操作主要通過高差閾值Δhi判斷[14]

(1)

其中:h0為最小高差閾值,s為研究區(qū)域的平均地形坡度,c為格網(wǎng)大小,hmax為最大高差閾值, 而wi為運(yùn)行到i次時(shí)濾波窗口的大小。當(dāng)前后格網(wǎng)內(nèi)的高程差小于高差閾值Δhi時(shí), 將當(dāng)前格網(wǎng)內(nèi)的點(diǎn)歸為地面點(diǎn)類格網(wǎng), 標(biāo)記為0; 而對(duì)于大于高差閾值的, 將該格網(wǎng)歸為非地面點(diǎn)類格網(wǎng), 標(biāo)記為1。通過遍歷所有點(diǎn)云數(shù)據(jù), 即可將地面點(diǎn)與非地面點(diǎn)分離開來。

2.2 基于區(qū)域增長與歐氏聚類相結(jié)合的點(diǎn)云分割

歐氏聚類算法是應(yīng)用于測(cè)繪領(lǐng)域的一種重要的點(diǎn)云分割算法。它利用KD樹搜索近鄰域點(diǎn), 計(jì)算鄰域點(diǎn)到該點(diǎn)的歐氏距離, 將在閾值范圍內(nèi)的點(diǎn)聚為一類, 通過反復(fù)迭代, 直到?jīng)]有新點(diǎn)加入為止。

(2)

其中:pi,qi∈Q,Q是一個(gè)點(diǎn)云集。具體的實(shí)現(xiàn)步驟為: ①對(duì)于空間中的某點(diǎn)p1, 利用KD樹搜索找到其近鄰點(diǎn), 計(jì)算這n個(gè)點(diǎn)到p1的距離, 將距離小于閾值的點(diǎn)p3、p4、p5…放到點(diǎn)集A中, 將點(diǎn)集A中的點(diǎn)從原始點(diǎn)云中移除; ②依次遍歷點(diǎn)集A中的所有點(diǎn), 重復(fù)步驟①, 更新點(diǎn)集A; ③當(dāng)A中無新點(diǎn)加入時(shí), 該類聚類完成; ④新建點(diǎn)集B, 對(duì)剩余點(diǎn)執(zhí)行步驟①操作; ⑤所有點(diǎn)都遍歷完成, 運(yùn)行結(jié)束。

通過歐氏聚類進(jìn)行分割的效果主要由設(shè)置的距離閾值決定: 當(dāng)設(shè)置的閾值過大時(shí), 會(huì)出現(xiàn)過分割現(xiàn)象, 建筑物與周邊較近的地物歸為一類,當(dāng)建筑物與樹木及地面周邊點(diǎn)粘連時(shí), 無法通過歐氏聚類分割(圖2a、2b); 反之, 則出現(xiàn)欠分割情形, 在圖2c、2d中, 由于閾值設(shè)置過小, 建筑物出現(xiàn)缺失現(xiàn)象。針對(duì)此問題并結(jié)合上述分析, 本文提出基于區(qū)域增長與歐氏聚類相結(jié)合的點(diǎn)云分割算法。

圖2 欠分割及過分割現(xiàn)象

區(qū)域增長算法是將點(diǎn)云的法向量進(jìn)行比較, 將滿足平滑約束條件的相鄰點(diǎn)歸為一類, 其工作原理為先利用KD樹對(duì)散亂點(diǎn)建立拓?fù)潢P(guān)系, 計(jì)算各點(diǎn)平均曲率值, 再將各點(diǎn)按曲率值從小到大進(jìn)行排序, 將曲率值最小的點(diǎn)設(shè)置為種子節(jié)點(diǎn), 由于曲率最小的點(diǎn)位于平坦區(qū)域, 從該區(qū)域進(jìn)行增長可以減少區(qū)域的總數(shù), 避免重復(fù)分割, 最后將滿足條件的不同地物聚為一類。

(3)

式中:Ki為曲面中點(diǎn)的平均曲率近鄰點(diǎn)；n為法向量;P為該點(diǎn)周圍無限小的區(qū)域; diam(P)為這個(gè)區(qū)域的直徑;為該點(diǎn)的梯度算子。該算法具體的步驟如下: ①將散亂點(diǎn)云按照曲率值進(jìn)行排序, 找到曲率最小的點(diǎn), 將其設(shè)置成種子節(jié)點(diǎn)，②計(jì)算該種子點(diǎn)與鄰近點(diǎn)的法線差值, 若小于設(shè)置的閾值, 且曲率值也小于設(shè)定的閾值, 則將該點(diǎn)添加到種子點(diǎn)集, 即歸為一類; ③將通過兩次檢驗(yàn)的點(diǎn)從原始點(diǎn)云中去除; ④設(shè)置最小和最大聚類的點(diǎn)數(shù); ⑤重復(fù)步驟①～③, 算法會(huì)將滿足條件的不同地物聚類為一類, 并用不同的顏色加以區(qū)分; ⑥當(dāng)剩余的點(diǎn)云數(shù)小于最小聚類點(diǎn)數(shù)時(shí), 則停止工作。

在大規(guī)模場(chǎng)景點(diǎn)云中, 建筑物等地物點(diǎn)云的曲率變化明顯要小于樹木等桿狀地物, 由于區(qū)域增長算法是按法線差值與曲率變化進(jìn)行分割的, 可將曲率變化較大的樹木及桿狀地物通過該算法分割出來, 歸為一類, 將剩余地物點(diǎn)云劃分為另一類, 再分別對(duì)兩類點(diǎn)云執(zhí)行歐氏聚類分割, 最終完成對(duì)典型地物的分割, 其分割流程見圖3。

圖3 典型地物分割流程圖

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

采用RIEGL VZ-1000掃描儀獲取的該數(shù)據(jù)，包含建筑物、周邊低矮樹木、路燈、旗桿等典型地物, 地物種類較多, 場(chǎng)景較為復(fù)雜。實(shí)驗(yàn)軟硬件環(huán)境:Windows 10專業(yè)版、處理器為Intel(R)Core(TM)i5-4200M CPU @2.50 GHz、內(nèi)存為8 GB; 本研究的實(shí)驗(yàn)算法是在Microsoft Visual Studio 2017編譯環(huán)境下, 結(jié)合PCL方法庫編程實(shí)現(xiàn)。

通過預(yù)處理后, 采用漸近形態(tài)學(xué)濾波算法移除地面點(diǎn), 濾波窗口設(shè)置為20 m, 坡度設(shè)置為1.0, 最小高差閾值設(shè)置為0.5 m, 最大高差閾值設(shè)置為0.8 m, 得到的結(jié)果如圖4所示。濾波后, 大量地面點(diǎn)都得以去除, 建筑物底部點(diǎn)及階梯處點(diǎn)未因該操作而受影響, 較好地保證了地物點(diǎn)的完整性。

圖4 地面點(diǎn)去除

經(jīng)形態(tài)學(xué)濾波后, 再對(duì)非地面點(diǎn)采用區(qū)域增長算法進(jìn)行分割, 法線差值設(shè)置為0.05, 曲率閾值設(shè)置為1.0。從圖5a可以看出, 由于樹及桿狀物曲率變化較大, 而建筑物點(diǎn)云分布較為平整, 因此樹及桿狀物聚類后呈現(xiàn)紅色, 而其余地物顏色與該類地物顏色不同；基于此，可將樹木及桿狀物點(diǎn)云簇分為一類(圖5b)；其余地物歸為另一類(圖5c)。接著，分別對(duì)上述兩類點(diǎn)云簇進(jìn)行歐氏聚類分割: 將典型地物個(gè)體單獨(dú)分離出來, 樹及桿狀物距離閾值設(shè)置成1.0 m, 最小聚類簇點(diǎn)數(shù)為50, 最大聚類簇點(diǎn)數(shù)為300 000, 最終生成的聚類簇?cái)?shù)為53個(gè)(圖6)；而對(duì)于建筑物等點(diǎn)云簇, 其距離閾值設(shè)置為2.0 m, 最小聚類簇點(diǎn)云為50, 最大聚類簇點(diǎn)數(shù)為300 000, 最后生成的聚類簇?cái)?shù)為7個(gè)(圖7)。

圖5 區(qū)域增長分割

圖6 樹及桿狀物歐氏聚類分割結(jié)果

圖7 建筑地物歐氏聚類結(jié)果

從圖8及圖9可以看出, 通過本文算法, 可將建筑物、樹、旗桿等典型地物單獨(dú)分割開來, 較好地解決了建筑物與地面、建筑物與樹桿、樹桿與地面等不同地物點(diǎn)云間因粘連問題而出現(xiàn)的過分割現(xiàn)象, 且本文將不同地物歸為兩類, 即樹及桿狀物點(diǎn)云簇和建筑物點(diǎn)云簇, 通過分別對(duì)其進(jìn)行歐氏聚類, 較好地解決了雙方對(duì)閾值設(shè)置的沖突, 避免出現(xiàn)地物欠分割現(xiàn)象。

圖8 樹及桿狀物分割

圖9 建筑物分割

4 結(jié) 論

以TLS獲取的大規(guī)模場(chǎng)景點(diǎn)云數(shù)據(jù)作為研究對(duì)象, 通過統(tǒng)計(jì)濾波與八叉樹精簡對(duì)源數(shù)據(jù)作預(yù)處理, 降低了噪聲點(diǎn)的干擾, 提高了點(diǎn)云分割的運(yùn)行效率。由于受地物形態(tài)特征、空間位置、地面點(diǎn)等因素影響, 在對(duì)典型地物進(jìn)行分割時(shí)，常會(huì)出現(xiàn)過分割或欠分割現(xiàn)象, 因此，本文提出基于區(qū)域增長與歐氏聚類相結(jié)合的分割算法, 首先通過形態(tài)學(xué)濾波去除地面點(diǎn), 然后通過區(qū)域增長算法, 將典型地物分為樹桿類及建筑物兩大類, 最后分別對(duì)這兩類地物通過歐氏聚類完成獨(dú)立個(gè)體的有效分割。本文方法有效解決了點(diǎn)云分割中易出現(xiàn)的過分割或欠分割問題, 在對(duì)樹、桿狀物及建筑物等典型地物的點(diǎn)云分割中有較強(qiáng)的實(shí)用性。然而, 本文算法的閾值需要人工調(diào)試, 如何自動(dòng)獲取合適的閾值, 將是下階段的研究重點(diǎn), 雖然本文事先對(duì)數(shù)據(jù)作了預(yù)處理, 提高了運(yùn)行效率, 但若數(shù)據(jù)量很大時(shí), 其運(yùn)行時(shí)間將成倍增長, 如何改進(jìn)算法, 加入并行策略, 減少分割時(shí)間, 將是下一步的工作重點(diǎn)。