基于最小二乘擬合的三維激光掃描點(diǎn)云濾波

嚴(yán)劍鋒，鄧喀中，邢正全

(1.國土環(huán)境與災(zāi)害監(jiān)測國家測繪地理信息局重點(diǎn)實(shí)驗(yàn)室，江蘇徐州221116;2.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇徐州221116)

一、引 言

地面三維激光掃描儀可快速獲取掃描區(qū)域的海量三維坐標(biāo)數(shù)據(jù)，在這些數(shù)據(jù)中，有進(jìn)行成圖、建模和分析等所需要的數(shù)據(jù)，也不可避免地存在一些影響測量結(jié)果的數(shù)據(jù)［1］。在地形掃描測量時(shí)，樹木植被的遮擋及長距離和大入射角掃描引起掃描點(diǎn)異常的可能性較大［2-3］。點(diǎn)云數(shù)據(jù)濾波需要去除這些噪聲點(diǎn)，只有將激光掃描點(diǎn)云數(shù)據(jù)進(jìn)行濾波后，才能有效地利用掃描數(shù)據(jù)，為最終生成DEM服務(wù)。目前，點(diǎn)云數(shù)據(jù)濾波方法主要包括基于地形坡度的方法、基于形態(tài)學(xué)開運(yùn)算方法、基于分類的算法、擬合法等［4-6］。

其中，擬合法以其原理簡單、去噪效果好的優(yōu)點(diǎn)而在地形掃描點(diǎn)云濾波方面得到廣泛應(yīng)用。傳統(tǒng)的移動(dòng)窗口曲面擬合濾波直接利用一個(gè)大尺度的移動(dòng)窗口通過尋找地面最低點(diǎn)擬合一個(gè)粗略地形模型，然后將高差超過給定閾值的點(diǎn)過濾掉［7］。這種方法雖然簡單，但過分依賴移動(dòng)窗口的大小。本文采用基于最小二乘的擬合濾波法，將二次曲面與多面函數(shù)結(jié)合。首先使用二次曲面擬合濾波，避免了直接進(jìn)行多面函數(shù)擬合錯(cuò)誤地選取噪聲點(diǎn)作為節(jié)點(diǎn)而對結(jié)果產(chǎn)生較大的影響;其次在二次曲面擬合去除大的噪聲點(diǎn)基礎(chǔ)上，再用多面函數(shù)擬合可以更好地逼近細(xì)節(jié)部分，選取準(zhǔn)確的地形點(diǎn)，進(jìn)而獲取準(zhǔn)確的DEM。

二、最小二乘擬合法去噪

二次曲面擬合可以總體上表達(dá)地形趨勢，而多面函數(shù)擬合的基本思想是：任何一個(gè)不規(guī)則的復(fù)雜曲面均可由一系列規(guī)則的數(shù)學(xué)表面總和以任意精度逼近［8-9］。因此，多面函數(shù)可以更好地顯示地形細(xì)節(jié)部分，但同時(shí)其擬合精度受節(jié)點(diǎn)的影響較大，這就需要保證擬合點(diǎn)的準(zhǔn)確性。本文將兩種擬合方法結(jié)合，首先用二次曲面擬合獲取大致趨勢面并去掉較大誤差點(diǎn)，為多面函數(shù)擬合提供精度較高的可選擬合點(diǎn);然后利用多面函數(shù)擬合濾波去除近地形表面的噪聲點(diǎn)。

最小二乘曲面擬合后需要設(shè)定閾值去除噪聲點(diǎn)，閾值的設(shè)定可根據(jù)濾波區(qū)域的地形情況及選用何種擬合方法確定。初次使用二次曲面擬合時(shí)，擬合面并不十分精確，只需去掉大的誤差點(diǎn)，閾值適當(dāng)取大，以免將非異常點(diǎn)去掉;第二次使用多面函數(shù)擬合去噪時(shí)，縮小閾值，去掉誤差較小的噪聲點(diǎn)，以獲取精確的點(diǎn)云。

1.二次曲面擬合［10］

對于高程變化平緩，即在沒有突變的范圍內(nèi)，可用曲面進(jìn)行逼近代替該區(qū)域。二次曲面擬合模型為

式中，(xk，yk，zk)為點(diǎn)的三維坐標(biāo);α0，α1，…，α5稱為二次曲面擬合系數(shù)。

當(dāng)擬合點(diǎn)數(shù)m大于必要測點(diǎn)數(shù)(二次曲面為6)時(shí)，有如下誤差公式

用最小二乘法可求得 α0，α1，…，α5的值，從而得到二次曲面表達(dá)式。

2.多面函數(shù)擬合［9，11］

設(shè)在測區(qū)內(nèi)有m個(gè)已測點(diǎn)S(x，y)，或記為數(shù)據(jù)點(diǎn)(x，y，s)，si為點(diǎn)(xi，yi)上的觀測量。用 n 個(gè)核函數(shù)的總和去逼近函數(shù)S(x，y)，即

式中，n為所取的已測點(diǎn)數(shù)，即節(jié)點(diǎn)數(shù)，n≤m;θ(x，y，xj，yj)為所取的核函數(shù);α為待定系數(shù)。

核函數(shù) θ(x，y，xj，yj)有多種類型可以選用，在測量數(shù)據(jù)處理中，一種常用的模型為具有對稱性的距離型倒雙曲面模型

式中，δ2為光滑因子，根據(jù)節(jié)點(diǎn)與擬合點(diǎn)的距離選用。

對于 m 個(gè)已知點(diǎn) Si(xi，yi)，i=1，2，…，m，方程矩陣形式為

其中，xjn、yjn表示所選的節(jié)點(diǎn)坐標(biāo)。由最小二乘法可求得系數(shù)α的值。

三、方法原理與實(shí)現(xiàn)

1.點(diǎn)云分割

掃描得到的點(diǎn)云數(shù)據(jù)量龐大，不論用何種方法進(jìn)行擬合去噪，都必須先選取一定數(shù)量的地表面點(diǎn)作為擬合點(diǎn)。將點(diǎn)云分成若干單元格，認(rèn)為每個(gè)單元格中的最低點(diǎn)為地面點(diǎn)，用最低點(diǎn)參與擬合。

按一定的長、寬、高將點(diǎn)云區(qū)域分割成m×n×k個(gè)基本小單元格。鄭德華［12］通過使用分割點(diǎn)云確定長方體單元格的方法進(jìn)行直接點(diǎn)云精簡。本文在實(shí)際的包含點(diǎn)云分割過程中，首先進(jìn)行二維網(wǎng)格劃分，再考慮第三維即z值的大小，以此得到需要的擬合點(diǎn)。具體的劃分步驟包括：

1)對于二維網(wǎng)格的劃分，為確定包含所有掃描點(diǎn)的長方形的大小，對于點(diǎn)云數(shù)據(jù)P(pi∈P)，搜索P中所有點(diǎn)的二維方向的最大值與最小值，即xmin、xmax、ymin、ymax。則點(diǎn)云數(shù)據(jù)二維邊界框的4個(gè)角點(diǎn)坐標(biāo)為：(xmin，ymin)、(xmin，ymax)、(xmax，ymin)、(xmax，ymax)。

2)將覆蓋整片點(diǎn)云的長方形區(qū)域沿x、y方向分割成m×n個(gè)小的基本長方形單元格。單元格的大小根據(jù)掃描區(qū)域及點(diǎn)云間隔來確定，在第一次粗略擬合去噪時(shí)一般單元格面積較大，網(wǎng)格邊長分別為gx、gy，即

3)確定每個(gè)單元格的邊界。對于每個(gè)單元格來說，需要給予相應(yīng)的編號，并確定該號的單元格的邊界，從而確定每個(gè)單元格內(nèi)的掃描點(diǎn)。邊界包括y方向的上界和下界、x方向的上界和下界，這樣就確定了每個(gè)長方形單元格所包含的所有點(diǎn)。

4)搜索每個(gè)單元格內(nèi)的高程最小點(diǎn)。掃描過程中由于遮擋或隨機(jī)噪聲的影響，點(diǎn)云分割后，不可避免地出現(xiàn)有些單元格為空的情況，因此人為去掉這些單元格是必要的。在剩余的單元格中逐個(gè)搜索高程最小點(diǎn)，得到一系列點(diǎn)，用于最小二乘擬合。

2.擬合去噪

用點(diǎn)云分割算法選取擬合點(diǎn)，利用上文中的方法進(jìn)行兩次擬合濾波。首先擬合二次曲面作為趨勢面，設(shè)定較大閾值，將大噪聲點(diǎn)去除;其次在初次濾波的基礎(chǔ)上再選取擬合點(diǎn)，并從中選取部分節(jié)點(diǎn)，進(jìn)行多面函數(shù)擬合，此時(shí)該多面函數(shù)可以很好地逼近地形表面，并精確地表達(dá)某些細(xì)節(jié)部分，去除近地表面噪聲點(diǎn)。

四、濾波試驗(yàn)和結(jié)果分析

本文掃描對象為一開采沉陷區(qū)域，大小約為85 m×14 m，分別在區(qū)域兩邊設(shè)站掃描，中間有部分重疊。此區(qū)域產(chǎn)生噪聲點(diǎn)的原因主要包括：①整片區(qū)域覆蓋了大面積雜草和植被，同時(shí)有些人為堆積的土坡，如圖1所示;②由于是在兩邊設(shè)站對中間區(qū)域掃描，因此中間部分得到的大入射角和長距離掃描點(diǎn)易成為噪聲點(diǎn);③一些隨機(jī)噪聲的影響。掃描共得到455 278個(gè)點(diǎn)。

圖1 原始點(diǎn)云

由原始掃描點(diǎn)云得到該區(qū)域的DEM，如圖2所示。由于噪聲點(diǎn)分布在整個(gè)區(qū)域，因此生成的DEM效果很差，區(qū)域高低沒有趨勢性，與實(shí)際不符。通過該原始點(diǎn)云數(shù)據(jù)無法得到準(zhǔn)確的DEM。

圖2 原始點(diǎn)云生成的DEM

進(jìn)行初次擬合濾波，通過點(diǎn)云單元格分割選取用于二次曲面擬合的點(diǎn)，本文實(shí)際選中的可用于擬合的點(diǎn)為1084個(gè)，得到如圖3所示的二次曲面。將二次曲面與剔除掉的噪聲點(diǎn)疊加(如圖4所示)，可以看出部分大的噪聲點(diǎn)已被去除。圖4中黑色長方形框中標(biāo)注的掃描時(shí)的標(biāo)靶，通過初次去噪已經(jīng)被部分剔除。

圖3 擬合的二次曲面(單位：m)

圖4 二次曲面與剔除的噪聲點(diǎn)疊加圖(單位：m)

從初次濾波后的點(diǎn)云中再選取擬合點(diǎn)，并從擬合點(diǎn)中選取多面函數(shù)擬合節(jié)點(diǎn)，進(jìn)行試驗(yàn)比較。由于擬合點(diǎn)較多，節(jié)點(diǎn)的分布對擬合的結(jié)果影響并不大，因此隨機(jī)選取在整個(gè)區(qū)域上均勻分布的節(jié)點(diǎn)。多面函數(shù)核函數(shù)選用距離型倒雙曲面模型，進(jìn)行多面函數(shù)擬合濾波。經(jīng)過兩次去噪后，剩余的點(diǎn)數(shù)為291 096，濾波后的點(diǎn)云側(cè)視如圖5所示。和原始點(diǎn)云相比，植被覆蓋及其他大量噪聲點(diǎn)已經(jīng)被剔除，濾波效果較為理想。

圖5 濾波后的點(diǎn)云側(cè)視圖

由濾波后的點(diǎn)云生成該區(qū)域DEM，如圖6所示。和圖2相比，該DEM較好地表達(dá)了實(shí)際地形趨勢。對濾波后準(zhǔn)確的DEM進(jìn)行三維渲染，結(jié)果如圖7所示。

圖6 濾波后點(diǎn)云生成的DEM

圖7 DEM三維渲染圖

對地形掃描的最終目的是建立精確的DEM及三維地形圖。由上述的實(shí)例分析可見，本文的濾波方法針對不大的區(qū)域點(diǎn)云濾波很有效，生成的DEM及三維渲染圖可以很好地反映該沉陷區(qū)的地形變化。

五、結(jié) 論

1)對于掃描得到的點(diǎn)云，考慮將誤差分為兩部分：較大誤差噪聲點(diǎn)和近地形表面噪聲點(diǎn)。前者包含了少量粗差;后者主要是由于隨機(jī)誤差和低矮植被遮擋引起的。這樣更有針對性地濾波可以取得較為理想的效果，進(jìn)而生成準(zhǔn)確的DEM。

2)提出點(diǎn)云分割算法，根據(jù)區(qū)域特點(diǎn)將點(diǎn)云區(qū)域首先進(jìn)行二維劃分，再考慮第三維，單元格形狀和大小的選取則視掃描區(qū)域范圍和地形而定，這可為兩次擬合提供滿足要求的擬合點(diǎn)。同時(shí)，由于區(qū)域長度較大，遠(yuǎn)距離處掃描點(diǎn)稀疏并且精度較低，若選取這些點(diǎn)作為擬合點(diǎn)，擬合出的曲面誤差較大，在點(diǎn)云分割時(shí)可將這些點(diǎn)去除，用附近的較為準(zhǔn)確點(diǎn)擬合的地形表面可以更好地表示該處地形。

3)傳統(tǒng)的移動(dòng)窗口最小二乘濾波只能適用于很小區(qū)域范圍去噪，同時(shí)過分依賴所取窗口大小。本文在此基礎(chǔ)上將兩種擬合方法結(jié)合，利用二次曲面和多面函數(shù)的不同特點(diǎn)對掃描區(qū)域進(jìn)行了擬合去噪，取得了較為理想的效果，為進(jìn)一步兩種或多種方法的組合濾波提供了借鑒。

［1］李亮，吳侃，劉虎，等.地面三維激光掃描地形測量數(shù)據(jù)粗差剔除算法及實(shí)現(xiàn)［J］.測繪科學(xué)，2010，35(3)：187-189.

［2］SOUDARISSANANE S，LINDENBERGH R，MENENTI M，et al.Scanning Geometry：Influenceing Factor on the Quality of Terrestrial Laser Scanning Points［J］.ISPRS Journal of Photogrammetry and Remote Sensing，2011，66(4)：389-399.

［3］張毅，閆利，楊紅，等.地面激光掃描球形標(biāo)靶的球心誤差研究［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2012，37(5)：598-601.

［4］VOSSELMAN G.Slope Based Filtering of Laser Altimetry Data［J］.IAPRS，2000(33)：935-942.

［5］ZHANG Keqi，CHEN Shu ching，WHITMAN D，et al.A Progressive Morphological Filter for Removing Nongroud Measurements from Airborne LIDAR Data［J］.IEEE Transactions on Geoscience and Remote Sensing，2003，41(4)：872-882.

［6］張皓，張永生，劉軍，等.一種基于平面擬合的LIDAR點(diǎn)云濾波方法［J］.測繪科學(xué)，2009，34(4)：141-143.

［7］PETZOLD B，REISSP，STOSSEL W.Laser Scaning-surveying and Mapping Agencies are Using a New Technique for the Deviation of Digital Terrain Models［J］.ISPRS Journal of Photogrammetry and Remote Sensing，1999，54(2-3)：95-104.

［8］呂言.數(shù)字地面模型中多面函數(shù)內(nèi)插法的研究［J］.武漢測繪學(xué)院學(xué)報(bào)，1981(2)：14-28.

［9］王新洲，陶本藻，邱衛(wèi)寧，等.高等測量平差［M］.北京：測繪出版社，2006.

［10］錢錦鋒.逆向工程中的點(diǎn)云處理［D］.杭州：浙江大學(xué)，2005.

［11］張菊清，劉平芝.抗差趨勢面與正交多面函數(shù)結(jié)合擬合 DEM 數(shù)據(jù)［J］.測繪學(xué)報(bào)，2008，37(4)：526-530.

［12］鄭德華.點(diǎn)云數(shù)據(jù)直接縮減方法及縮減效果研究［J］.測繪工程，2006，15(4)：27-30.