第python點云地面點濾波目錄1.引言1.1什么是地面點濾波？1.2地面點濾波的相關方法2.PTD具體介紹2.1去除孤立點2.2參數說明2.3選擇種子點并構建TIN模型2.4迭代加密TIN模型本篇博客參考：

1）DEMgenerationfromlaserscannerdatausingadaptiveTINmodels

2）FilteringairborneLiDARdatabyembeddingsmoothness-constrainedsegmentationinprogressiveTINdensification

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

1.引言

1.1什么是地面點濾波？

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

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

1.2地面點濾波的相關方法

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

slope-based

block-minimum

surface-based

clustering/segmentationalgorithm

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

Morphology-basedfilters

Iterative-interpolation-basedfilters

Progressive-densification-basedfilters

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

2.PTD具體介紹

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

2.1去除孤立點

如引言中1.1所述，去除孤立點類似于圖像中的去噪操作。outliers是測量數據集中的那些遠高于/低于地表的點，這中情況常常會導致濾波算法出現錯誤（例如，PTD算法中假設格網中的最低點為地面點，從而導致錯誤）。可以通過下述三個簡單步驟來去除孤立點，當然除此之外還有其他很多更加優(yōu)秀的方法。

對所有數據的高程進行統(tǒng)計，建立一個高程分布的直方圖，通過觀察高程分布來確定高程閾值，從而消除分布中最低和最高的小“尾巴”。

通過每個點與周圍點之間的最小高程差來進行搜尋仍然存在的孤立點（這里使用一個2D的KD樹來進行組織查詢每個點的近鄰點）。

手工校正孤立點自動去除過程中所產生的錯誤。

2.2參數說明

在PTD算法中有以下6個參數進行預先設置：

1）最大建筑尺寸m：m是一個長度閾值，此閾值被用來定義格網的大小，隨后算法可以處理建筑物尺寸小于此閾值的建筑物。

2）最大地形角度t：t是一個坡度閾值，決定了通過什么方式（是否進行設置鏡像點）去判斷未分類點的類別。如果未分類點所在三角面的坡度大于t則應該通過一個鏡像點來進行判斷，反之則直接判斷。（后續(xù)在2.4中也會再進行詳細介紹）。

3）最大角度θ：θ是三角面與待分類點和最近的三角網頂點之間連線之間的最大角度。如果一個未分類點對應的角度大于θ則被標記為地物點，否則設置為地面點。

4）最大距離d：d是當前迭代中從待判斷點到三角面之間的最大距離，類似的，如果一個未分類點對應的最大距離大于d則被標記為地物點，否則設置為地面點。

5）最小邊長l：l是構建TIN模型中所有三角形最長邊（平面投影）的最小閾值。當三角形中的所有邊都小于l時，則停止在三角網中加入地面點（注意l是在平面中計算的）。因此，此參數可以避免引起地面模型中過高的點密度，以及降低內存的使用。

6）最大邊長l'：*l'是構建TIN模型中所有三角形最短邊（平面投影）的最大閾值，當三角形中的所有邊都小于l'*被用于停止處理處理三角形。因此，此參數用于稀疏地面點，以及降低內存的使用。

2.3選擇種子點并構建TIN模型

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

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

根據整個數據集劃分為格網之后，每個網格中的最低點被設置為“種子點”（初始地面點）。除此之外，“boundingbox”的四個角點也被設置為“種子點”（其高程值等同于距離最近的種子點高程），如下圖所示：

注：把“boundingbox”的四個角點也被設置為“種子點”是為了保證所有點都處于TIN模型內部。

隨后，根據選擇好的種子點來構建初始TIN模型來表示初始地表，剩余的點被默認標記為地物。

2.4迭代加密TIN模型

在每次迭代過程中通過預先設置的閾值參數，來對“潛在點(potentialpoint)”進行逐點判斷。詳細步驟如下所述：

1）確定潛在點(potentialpoint)的所在位置PPotential(xp,yp,zp)，找到其所在的三角形Ttriangle，即PPotential在三角形內部或邊緣上或者在頂點上。

2）計算三角平面的坡度Striangle，如果Striangle小于預先設置的最大地形角度t，則進行第3)步，若大于則進行第4)步。

3）如下圖所示，計算的兩個參數，包括：三角面Ttriangle與待分類點PPotential和最近的三角網頂點之間連線之間的角度，表示為Aangle，以及待分類點PPotential和三角面Ttriangle之間的距離，表示為Ddistance。如果PPotential所對應的上述兩個參數Aangle小于最大角度θ，Ddistance小于最大距離d則認為待分類點PPotential是地面點，否則認為是地物點。隨后，進行下一個點的判斷。

4）如果三角平面的坡度Striangle大于預先設置的最大地形角度t則需要設置PPotential的鏡像點。先找到PPotential所在三角形

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

其中(xmirror,ymirror,zmirror)是所求鏡像點的三維坐標。

接著對鏡像點使用步驟3）的方式來進行計算Aangle與Ddistance兩個參數進行判斷，來決定PPotential的類別。

5）在結束每次迭代之后，新檢測出的地面點通過下述步驟加入TIN模型中。

確定Pground(xg,yg,zg)的坐標，找到其所在的三角形T'triangle。

計算三角形T'triangle的所有邊在水平投影中的長度，如果任意邊的長度都大于