python 點(diǎn)云地面點(diǎn)濾波-progressive TIN densification(PTD)算法介紹

本篇博客參考：
1）DEM generation from laser scanner data using adaptive TIN models
2）Filtering airborne LiDAR data by embedding smoothness-constrained segmentation in progressive TIN densification

文章名中有超鏈接，若不方便下載，則可以在此：資源鏈接進(jìn)行下載。

1.引言

1.1什么是地面點(diǎn)濾波？

機(jī)載激光雷達(dá)(airborne light detection and ranging)/機(jī)載激光掃描(ALS, airborne laser scanning)在過去20多年的時(shí)間里迅速發(fā)展，其相對于傳統(tǒng)攝影測量影像及INSAR(干涉合成孔徑雷達(dá))可以直接記錄從地物或地表返回的密集、離散、細(xì)節(jié)豐富、精確的三維點(diǎn)云。如何對這些不規(guī)則點(diǎn)云進(jìn)行處理應(yīng)用是我們要解決的問題，其中一個(gè)重要的處理步驟就是：地面點(diǎn)濾波(ground filterring)，簡而言之就是“在無序、不規(guī)則的三維離散點(diǎn)云中找到哪些是由地表返回的，哪些是由地物返回的?！?/p>

注：關(guān)于地面點(diǎn)濾波的概念我們要與孤立點(diǎn)(outlier)濾波區(qū)分開，孤立點(diǎn)濾波可以理解為圖像中的去噪，去除數(shù)據(jù)測量過程中受到飛鳥、多路徑效應(yīng)所產(chǎn)生的遠(yuǎn)低于/高于其他數(shù)據(jù)的點(diǎn)。

1.2地面點(diǎn)濾波的相關(guān)方法

眾多學(xué)者已經(jīng)提出了各種類型的濾波算法來從ALS三維點(diǎn)云中自動提取地面點(diǎn)，根據(jù)濾波器的概念可以分為以下四類：

• slope-based
• block-minimum
• surface-based
• clustering/segmentation algorithm

基于面的地面點(diǎn)濾波方法的核心步驟是創(chuàng)建一個(gè)最接近裸露地表的表面，其使用了更多的context(上下文，環(huán)境、背景)信息，所以一般可以取得比其他濾波方法更好的濾波效果。此外，根據(jù)創(chuàng)建表面的方法又可以把surface-based類型濾波器分為以下三個(gè)子類：

• Morphology-based filters
• Iterative-interpolation-based filters
• Progressive-densification-based filters

其中Morphology-based filters使用不同尺寸大小的窗口形態(tài)學(xué)操作(opening/geodesic)來去除不同大小的地物，但是這類方法需要假設(shè)地形的坡度是一個(gè)常數(shù)。以及面臨著一個(gè)巨大的挑戰(zhàn)是在窗口大小變化的時(shí)候怎么保持地形特征不變；Iterative-interpolation-based filters是通過整個(gè)點(diǎn)云數(shù)據(jù)集來逐步接近地表，首先通過一個(gè)粗糙的表面來計(jì)算點(diǎn)云到表面的殘差，通常來說地物點(diǎn)會有正的差值，地面點(diǎn)會有負(fù)的差值。這種方法最大的挑戰(zhàn)是當(dāng)精度一定的條件下怎么來提高算法的效率；Progressive-densification-based filters與上一種方法類似，也是漸進(jìn)地把每個(gè)點(diǎn)逐步的分類為地面點(diǎn)，然而此方法不需要進(jìn)行插值。本篇博客中所介紹的PTD(progressive TIN densification)就是屬于surface-based類型中的一種Progressive-densification-based filters。

2.PTD具體介紹

PTD是Axelsson在2000年左右提出的一種經(jīng)典濾波方法，在工程應(yīng)用（TerraScan）及科學(xué)社區(qū)中得到了廣泛的應(yīng)用。主要可以通過以下五個(gè)步驟進(jìn)行實(shí)現(xiàn)：

2.1去除孤立點(diǎn)

如引言中1.1所述，去除孤立點(diǎn)類似于圖像中的去噪操作。outliers是測量數(shù)據(jù)集中的那些遠(yuǎn)高于/低于地表的點(diǎn)，這中情況常常會導(dǎo)致濾波算法出現(xiàn)錯(cuò)誤（例如，PTD算法中假設(shè)格網(wǎng)中的最低點(diǎn)為地面點(diǎn)，從而導(dǎo)致錯(cuò)誤）?？梢酝ㄟ^下述三個(gè)簡單步驟來去除孤立點(diǎn)，當(dāng)然除此之外還有其他很多更加優(yōu)秀的方法。

• 對所有數(shù)據(jù)的高程進(jìn)行統(tǒng)計(jì)，建立一個(gè)高程分布的直方圖，通過觀察高程分布來確定高程閾值，從而消除分布中最低和最高的小“尾巴”。
• 通過每個(gè)點(diǎn)與周圍點(diǎn)之間的最小高程差來進(jìn)行搜尋仍然存在的孤立點(diǎn)（這里使用一個(gè)2D的KD樹來進(jìn)行組織查詢每個(gè)點(diǎn)的近鄰點(diǎn)）。
• 手工校正孤立點(diǎn)自動去除過程中所產(chǎn)生的錯(cuò)誤。

 2.2參數(shù)說明

在PTD算法中有以下6個(gè)參數(shù)進(jìn)行預(yù)先設(shè)置：

1）最大建筑尺寸m：m是一個(gè)長度閾值，此閾值被用來定義格網(wǎng)的大小，隨后算法可以處理建筑物尺寸小于此閾值的建筑物。
2）最大地形角度t：t是一個(gè)坡度閾值，決定了通過什么方式（是否進(jìn)行設(shè)置鏡像點(diǎn)）去判斷未分類點(diǎn)的類別。如果未分類點(diǎn)所在三角面的坡度大于t則應(yīng)該通過一個(gè)鏡像點(diǎn)來進(jìn)行判斷，反之則直接判斷。（后續(xù)在2.4中也會再進(jìn)行詳細(xì)介紹）。
3）最大角度θ：θ是三角面與待分類點(diǎn)和最近的三角網(wǎng)頂點(diǎn)之間連線之間的最大角度。如果一個(gè)未分類點(diǎn)對應(yīng)的角度大于θ則被標(biāo)記為地物點(diǎn)，否則設(shè)置為地面點(diǎn)。
4）最大距離d：d是當(dāng)前迭代中從待判斷點(diǎn)到三角面之間的最大距離，類似的，如果一個(gè)未分類點(diǎn)對應(yīng)的最大距離大于d則被標(biāo)記為地物點(diǎn)，否則設(shè)置為地面點(diǎn)。
5）最小邊長l：l是構(gòu)建TIN模型中所有三角形最長邊（平面投影）的最小閾值。當(dāng)三角形中的所有邊都小于l時(shí)，則停止在三角網(wǎng)中加入地面點(diǎn)（注意l是在平面中計(jì)算的）。因此，此參數(shù)可以避免引起地面模型中過高的點(diǎn)密度，以及降低內(nèi)存的使用。
6）最大邊長l'：*l'是構(gòu)建TIN模型中所有三角形最短邊（平面投影）的最大閾值，當(dāng)三角形中的所有邊都小于l'*被用于停止處理處理三角形。因此，此參數(shù)用于稀疏地面點(diǎn)，以及降低內(nèi)存的使用。

2.3選擇種子點(diǎn)并構(gòu)建TIN模型

對給定的點(diǎn)云數(shù)據(jù)集定義一個(gè)特定的“bounding box”并固定左上角坐標(biāo)(xtopleft, ytopleft)右下角坐標(biāo)(xbottomright, ybottomright)，寬度w，以及高度h。然后通過上述定義的變量通過以下公式把整個(gè)數(shù)據(jù)區(qū)域劃分成nRow行，nColumn列，尺寸大小為m的格網(wǎng)。

其中m，為最大建筑尺寸，ceil(x)函數(shù)代表向上取整，即找到不小于x的最小整數(shù)。

根據(jù)整個(gè)數(shù)據(jù)集劃分為格網(wǎng)之后，每個(gè)網(wǎng)格中的最低點(diǎn)被設(shè)置為“種子點(diǎn)”（初始地面點(diǎn)）。除此之外，“bounding box”的四個(gè)角點(diǎn)也被設(shè)置為“種子點(diǎn)”（其高程值等同于距離最近的種子點(diǎn)高程），如下圖所示：

注：把“bounding box”的四個(gè)角點(diǎn)也被設(shè)置為“種子點(diǎn)”是為了保證所有點(diǎn)都處于TIN模型內(nèi)部。

隨后，根據(jù)選擇好的種子點(diǎn)來構(gòu)建初始TIN模型來表示初始地表，剩余的點(diǎn)被默認(rèn)標(biāo)記為地物。

2.4迭代加密TIN模型

在每次迭代過程中通過預(yù)先設(shè)置的閾值參數(shù)，來對“潛在點(diǎn)(potential point)”進(jìn)行逐點(diǎn)判斷。詳細(xì)步驟如下所述：
1）確定潛在點(diǎn)(potential point)的所在位置PPotential(xp,yp,zp)，找到其所在的三角形Ttriangle，即PPotential在三角形內(nèi)部或邊緣上或者在頂點(diǎn)上。
2）計(jì)算三角平面的坡度Striangle，如果Striangle小于預(yù)先設(shè)置的最大地形角度t，則進(jìn)行第3)步，若大于則進(jìn)行第4)步。
3）如下圖所示，計(jì)算的兩個(gè)參數(shù)，包括：三角面Ttriangle與待分類點(diǎn)PPotential和最近的三角網(wǎng)頂點(diǎn)之間連線之間的角度，表示為Aangle，以及待分類點(diǎn)PPotential和三角面Ttriangle之間的距離，表示為Ddistance。如果PPotential所對應(yīng)的上述兩個(gè)參數(shù)Aangle小于最大角度θ，Ddistance小于最大距離d則認(rèn)為待分類點(diǎn)PPotential是地面點(diǎn)，否則認(rèn)為是地物點(diǎn)。隨后，進(jìn)行下一個(gè)點(diǎn)的判斷。

4）如果三角平面的坡度Striangle大于預(yù)先設(shè)置的最大地形角度t則需要設(shè)置PPotential的鏡像點(diǎn)。先找到PPotential所在三角形

Ttriangle中高程值最大的頂點(diǎn)Pvertex(xv,yv,zv)，然后通過下式計(jì)算PPotential的鏡像點(diǎn)：

其中(xmirror,ymirror,zmirror)是所求鏡像點(diǎn)的三維坐標(biāo)。

接著對鏡像點(diǎn)使用步驟3）的方式來進(jìn)行計(jì)算Aangle與Ddistance兩個(gè)參數(shù)進(jìn)行判斷，來決定PPotential的類別。
5）在結(jié)束每次迭代之后，新檢測出的地面點(diǎn)通過下述步驟加入TIN模型中。

• 確定Pground(xg,yg,zg)的坐標(biāo)，找到其所在的三角形T'triangle。
• 計(jì)算三角形T'triangle的所有邊在水平投影中的長度，如果任意邊的長度都大于l加入當(dāng)前地面點(diǎn)Pground到TIN模型中，并進(jìn)行刷新。否則，進(jìn)行判斷下一個(gè)新檢測的地面點(diǎn)。
• 重復(fù)上述迭代，直到不再有點(diǎn)被加入到地面點(diǎn)集中。

上述五個(gè)步驟就是PTD算法的主要實(shí)現(xiàn)方法了，PTD算法已經(jīng)被廣泛應(yīng)用在各種類型的景觀（地形）中，且取得不錯(cuò)的效果。但是需要注意的是，雖然在算法中采用了設(shè)置“鏡像點(diǎn)”的方式來避免出現(xiàn)cutting-off的問題（斷裂線分布區(qū)域），其還是對陡峭地形非常敏感。

到此這篇關(guān)于python 點(diǎn)云地面點(diǎn)濾波-progressive TIN densification(PTD)算法介紹的文章就介紹到這了,更多相關(guān)python PTD點(diǎn)云地面點(diǎn)濾波內(nèi)容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！