基于最小二乘擬合LiDAR點(diǎn)云濾波改進(jìn)算法

周長江+何永鐘

摘 要：針對(duì)傳統(tǒng)的使用單一擬合法濾波的不足，提出組合曲面擬合濾波方法。通過引入虛擬網(wǎng)格點(diǎn)方法，改進(jìn)數(shù)據(jù)預(yù)處理方法，在剔除掉粗差點(diǎn)云后，選用較為準(zhǔn)確的多面函數(shù)擬合，更好的表達(dá)地形起伏細(xì)節(jié)，逼近真實(shí)地形表面，進(jìn)行二次濾波。實(shí)例證明與ISPR公布的經(jīng)典濾波算法相比，該組合最小二乘擬合的點(diǎn)云濾波效果較好，可用于生成準(zhǔn)確的DEM。

關(guān)鍵詞：點(diǎn)云數(shù)據(jù) 多面函數(shù) 最小二乘 濾波

中圖分類號(hào)：P237 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1672-3791（2017）05（b）-0240-03

LiDAR點(diǎn)云除去地面點(diǎn)外，還含有樹木、建筑物、橋梁、車輛等信息，真實(shí)的DEM數(shù)字高程模型需要過濾掉非地面點(diǎn)。目前，點(diǎn)云過濾方法主要基于形態(tài)學(xué)、基于坡度、基于內(nèi)插、基于移動(dòng)曲面擬合等算法。在總結(jié)對(duì)比分析這些算法基礎(chǔ)上，該文通過將最小二乘與多面函數(shù)擬合兩種方式相結(jié)合，目的為嘗試獲取更加準(zhǔn)確的DEM。

1 最小二乘擬合法濾波

二次曲面擬合可以從總體上表示地形走勢，而多面函數(shù)能較好地體現(xiàn)地形細(xì)節(jié)，但由于其擬合精度受節(jié)點(diǎn)影響較大。該文通過將兩種方法結(jié)合，首先對(duì)點(diǎn)云進(jìn)行預(yù)處理，再將二次曲面擬合獲取大致趨勢面并去掉較大誤差點(diǎn)，為多面函數(shù)擬合提供精度較高的可選擬合點(diǎn)；然后利用多面函數(shù)擬合濾波去除非地面點(diǎn)。

1.1 二次曲面擬合

針對(duì)高程變化不大區(qū)域，可用二次曲面進(jìn)行模擬，擬合之?dāng)?shù)學(xué)模型為[2]：

2.2 多面函數(shù)擬合

設(shè)在測區(qū)域內(nèi)有m個(gè)已經(jīng)測點(diǎn)云S（x，y），標(biāo)記成數(shù)據(jù)點(diǎn)（x，y，s），Si為點(diǎn)（xi，yi）上的觀側(cè)值。用n個(gè)核函數(shù)的總和去逼近函數(shù)S（x，y）[3]，即：

其中，表達(dá)選擇地節(jié)點(diǎn)坐標(biāo)，運(yùn)用最小二乘法可以求出系數(shù)α的值。

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

2.1 點(diǎn)云數(shù)據(jù)預(yù)處理

該文通過對(duì)單一閾值法進(jìn)行改進(jìn)，目的是在數(shù)據(jù)預(yù)處理環(huán)節(jié)剔除低位誤差。數(shù)據(jù)處理過程包含兩步：（1）首先假定一個(gè)窗口區(qū)域，將該區(qū)域內(nèi)幾何中心點(diǎn)與其他的激光腳點(diǎn)進(jìn)行高程跳變檢測，對(duì)比高程的最小差值與所設(shè)置的閾值，當(dāng)中心點(diǎn)與其他點(diǎn)的最小差值大于閾值時(shí)，認(rèn)定為是粗差點(diǎn)；（2）第二步繼續(xù)對(duì)點(diǎn)云數(shù)據(jù)集進(jìn)行搜索，通過設(shè)定閾值作為判定條件，如果搜索到的點(diǎn)相對(duì)較集中并且符合判定條件，則認(rèn)定該點(diǎn)集內(nèi)的點(diǎn)為粗差點(diǎn)。

2.2 點(diǎn)云格網(wǎng)化分割

機(jī)載LiDAR點(diǎn)云數(shù)據(jù)量比較龐大，曲面擬合時(shí)，需要先選取一定量的地面點(diǎn)作為擬合點(diǎn)。首先將預(yù)處理后的點(diǎn)云進(jìn)行格網(wǎng)化，分隔成若干個(gè)單元格，用每個(gè)單元格中最點(diǎn)作為地面點(diǎn)參與擬合。具體的劃分步驟包括：

（1）對(duì)于點(diǎn)云數(shù)據(jù)，搜索Q中所有點(diǎn)的二維方向的最大值最小值，即、、、。

（2）單元格之劃分應(yīng)根據(jù)掃描區(qū)域點(diǎn)云的密度D來確定，單個(gè)正方形單元格的最低長度應(yīng)為 m。對(duì)LiDAR點(diǎn)云數(shù)據(jù)進(jìn)行格網(wǎng)化，沿x，y方向分割成mm個(gè)小的單元，網(wǎng)格邊長分別為，即：

（3）確定每個(gè)單元格的邊界。通過定義單元格編號(hào)來明確各單元格邊界及范圍內(nèi)需要掃描點(diǎn)。邊界涵蓋了x、y兩個(gè)方向的上下界。

（4）逐個(gè)搜索單元格內(nèi)高程值最小點(diǎn)。在對(duì)點(diǎn)云進(jìn)行網(wǎng)格分割之后，會(huì)不可避免地出現(xiàn)部分單元格為空值的情況，所以人為去掉這些單元格是十分必要的。然后，在剩余的單元格中逐個(gè)搜索高程值最小點(diǎn)，用于最小二乘擬合。

2.3 自適應(yīng)閾值確定

二次曲面擬合之后，原始地面點(diǎn)高程值與擬合高程值是相對(duì)比較接近的，對(duì)于非地面點(diǎn)而言仍具有較大的高差。該文選用一種局部自適應(yīng)閾值方法，閾值的選取采用每個(gè)單元網(wǎng)格內(nèi)高程最大值減去最低點(diǎn)后乘以一個(gè)系數(shù)。即：

Factoy的取值由格網(wǎng)內(nèi)地形起伏決定，即統(tǒng)計(jì)每個(gè)格網(wǎng)內(nèi)高程數(shù)據(jù)求算平均坡度。

2.4 曲面擬合濾波

采用點(diǎn)云格網(wǎng)化分割算法選取適量擬合點(diǎn)，基于前面流程進(jìn)行兩次擬合濾波。首先將擬合二次曲面作為趨勢面，采用自適應(yīng)閾值，剔除高差大的點(diǎn)； 然后在第一次濾波的基礎(chǔ)上再選取部分節(jié)點(diǎn)，進(jìn)行多面函數(shù)擬合，此時(shí)該多面函數(shù)可以精確地表達(dá)某些細(xì)節(jié)部分，濾除非地面點(diǎn)。

3 濾波試驗(yàn)與結(jié)果分析

該文的實(shí)驗(yàn)數(shù)據(jù)來自2003年ISPRS （國際攝影測量與遙感協(xié)會(huì)）第Ⅲ工作小組提供的比較實(shí)驗(yàn)數(shù)sample12，該樣例數(shù)據(jù)掃描區(qū)域內(nèi)地形相對(duì)比較復(fù)雜點(diǎn)，建筑物及樹木、植被地形要素較多。區(qū)域大小約為204 m×264 m，點(diǎn)云密度為0.67 points/，點(diǎn)云總數(shù)為52 119?；谠键c(diǎn)云數(shù)據(jù)通過內(nèi)插生成DEM，見圖1。因建筑物、樹木、植被等噪聲點(diǎn)分布于整個(gè)測區(qū)，生成的DEM精度較差，與實(shí)際地形不符。

利用直方圖對(duì)原始地點(diǎn)云高程值進(jìn)行數(shù)理統(tǒng)計(jì)（見圖2），從樣例數(shù)據(jù)點(diǎn)云高程直方圖上可以直觀看出：Sample12存在著較為明顯低位粗差點(diǎn)。因此采用改進(jìn)的單一閾值法先對(duì)樣例數(shù)據(jù)進(jìn)行預(yù)處理，剔除低位粗差點(diǎn)（見圖2）。參數(shù)設(shè)置為：窗口半徑設(shè)置為10 m、高程閾值設(shè)值為2 m、粗差點(diǎn)集最大允許包含數(shù)為6個(gè)，共有25個(gè)粗差點(diǎn)從樣例數(shù)據(jù)中剔除。

由于兩個(gè)激光點(diǎn)距離過遠(yuǎn)時(shí)，對(duì)擬合值的影響較小，且擬合點(diǎn)數(shù)過多時(shí)，一定程度上也會(huì)影響效率。因此，實(shí)際處理時(shí)也要先對(duì)點(diǎn)云坐標(biāo)進(jìn)行統(tǒng)計(jì)，進(jìn)而將點(diǎn)云集合分隔成10×10的格網(wǎng)，選擇格網(wǎng)中最低點(diǎn)參與初次擬合。該文實(shí)際用于擬合的點(diǎn)數(shù)為100，得到如圖3所示的二次曲面。通過得到的二次曲面，根據(jù)每個(gè)格網(wǎng)地形坡度設(shè)定相應(yīng)的閾值，濾除大部分的建筑物及植被等非地面點(diǎn)。

點(diǎn)云完成第一次擬合過濾后，從剩余點(diǎn)云集合中，選取多面函數(shù)擬合節(jié)點(diǎn)，進(jìn)行試驗(yàn)比較。由于擬合點(diǎn)數(shù)量相對(duì)較多，節(jié)點(diǎn)的位置分布對(duì)擬合的結(jié)果影響并不大，因此可以在區(qū)域內(nèi)隨機(jī)選取分布較為均勻的節(jié)點(diǎn)。多面函數(shù)核函數(shù)選用距離型倒雙曲面模型，進(jìn)行多面函數(shù)擬合濾波。經(jīng)過兩次去噪后，剩余的點(diǎn)數(shù)為26 144。由濾波后的點(diǎn)云生成該區(qū)域DEM，該文還采用了定量分析的方法進(jìn)行該濾波算法的評(píng)價(jià)。機(jī)載激光雷達(dá)數(shù)據(jù)濾波中，存在兩類誤差：第Ⅰ類誤差是地面點(diǎn)被誤分為非地面點(diǎn)的百分比，第Ⅱ類誤差是非地面點(diǎn)被誤分為地面點(diǎn)的百分比，第Ⅲ類誤差為總誤差，即被錯(cuò)分的點(diǎn)占整個(gè)數(shù)據(jù)的百分比。統(tǒng)計(jì)選用的sample12濾波后各類誤差，將各類誤差與ISPRS公布的8類經(jīng)典算法做了比較，見表1。

由表1分析可得，該文濾波算法與8類經(jīng)典算法相比：第Ⅰ類誤差顯著減小，第Ⅱ類誤差也不會(huì)明顯增大，同時(shí)總誤差也都相應(yīng)有所減小，濾波結(jié)果較好，生成的DEM也能很好地反映區(qū)域地形變化。

4 結(jié)論

（1）對(duì)于傳統(tǒng)的移動(dòng)曲面最小二乘濾波只適用于較小區(qū)域?yàn)V波，而且過分依賴所取曲面種子點(diǎn)，該文將兩種擬合方法相結(jié)合，利用二次曲面和多面函數(shù)的不同特點(diǎn)對(duì)掃描區(qū)域進(jìn)行了擬合去噪，取得了較為理想的效果。

（2） 改進(jìn)了單一閾值法，剔除了低位粗差點(diǎn)的影響；并改進(jìn)了點(diǎn)云格網(wǎng)化分割算法，根據(jù)區(qū)域特點(diǎn)將區(qū)域先進(jìn)行二維劃分，再考慮第三維，單元格形狀和大小的選取則視掃描區(qū)域范圍而定，這可為兩次擬合提供滿足要求的擬合點(diǎn)。

參考文獻(xiàn)

[1] 嚴(yán)劍鋒，鄧喀中，邢正全.基于最小二乘擬合的三維激光掃描點(diǎn)云濾波[J].測繪通報(bào)，2013（5）：43-46.

[2] 房華樂.面向?qū)ο蟮能囕dLiDAR點(diǎn)云濾波方法[J].測繪科學(xué)，2015（4）：92-96

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

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