点云地面滤波

形态学滤波主要包括腐蚀和膨胀以及二者相结合产生的开、闭操作。 腐蚀：即去除不必要的部分，简化物体的形状。举个例子：一棵大树，去掉树叶，只保留树干，用树干表示树木，即提取树木的“骨架”，保留主要信息，表示原来的物体。可以去掉一些噪点，简化点云。 膨胀：即在物体原有的基础上，增加物体的体积。举个例子：一棵大树，使得树叶增加茂密，树干更加粗壮，在视觉上显得树木细节更加丰富，依次表示原来的物体。可以修复一些空洞使得细节更加饱满。 开操作：先腐蚀后膨胀的操作称之为开操作。它具有消除细小物体，在纤细处分离物体和平滑较大物体边界的作用，使得物体的骨架更加突出。 闭操作：先膨胀后腐蚀的操作称之为闭操作。它具有填充物体内细小空洞，连接邻近物体和平滑边界的作用，同样可以使得物体骨架更加突出。

窗口大小对于形态学地面滤波至关重要。因此形态学地面滤波重点讨论如何确定最优的窗口大小。对于这个问题可以通过逐渐增大形态学滤波器的窗口尺寸来解决，这种方法被称为渐进式形态学滤波。下图说明这种方法的过程：

上图表示渐进式形态学滤波识别地形和建筑物量测的过程。这些点代表基于LIDAR采集的点云。第一个滤波高程面(虚线)是通过对原始点数据应用窗口大小为15 m 的开运算得到的。通过对第一个滤波曲面施加窗口大小为21 m的开运算得到第二个滤波高程曲面(实线)。 然而，在上图中的滤波过程往往会产生一个位于地形测量点云下方的表面，导致高处起伏地形顶部的测量点云被错误去除。即使在平坦的地面区域，过滤后的表面通常位于原始测量点云下方。因此，对于地形的大部分测量点云会被去除。这个问题可以通过引入基于地形、建筑物和树木的高度变化的高度差阈值来克服。下图说明了渐进式形态学滤波的主要过程：

在应用形态学滤波时，窗口大小和高差阈值的选取对取得良好的效果至关重要。对于窗口大小的选择，一个直观的选择是通过以下公式线性增加窗口大小 w k w_k wk​： w k = 2 k b + 1 (1) w_k=2 k b+1\tag{1} wk​=2kb+1(1) k k k:迭代次数， b b b:初始窗口大小。然而，对于具有大型非地面物体的区域，需要相当长的计算时间。 或者也可以采用下面这种方法： w k = 2 b k + 1 (2) w_k=2 b^k+1\tag{2} wk​=2bk+1(2)

高差阈值可根据研究区地形坡度确定。假设坡度不变，地形最大高差 d h max ⁡ ( t ) , k d h_{\max (t), k} dhmax(t),k​、窗口大小 w k w_k wk​与地形坡度 s s s之间存在关系： s = d h max ⁡ ( t ) , k ( w k − w k − 1 ) 2 (3) s=\frac{d h_{\max (t), k}}{\frac{\left(w_k-w_{k-1}\right)}{2}}\tag{3} s=2(wk​−wk−1​)​dhmax(t),k​​(3) 高差阈值 d h T , k d h_{T, k} dhT,k​: d h T , k = { d h 0 ,  if  w k ≤ 3 s ( w k − w k − 1 ) c + d h 0 ,  if  w k > 3 d h max ⁡ ,  if  d h T , k > d h max ⁡ (4) d h_{T, k}= \begin{cases}d h_0, & \text { if } w_k \leq 3 \\ s\left(w_k-w_{k-1}\right) c+d h_0, & \text { if } w_k>3 \\ d h_{\max }, & \text { if } d h_{T, k}>d h_{\max }\end{cases}\tag{4} dhT,k​=⎩ ⎨ ⎧​dh0​,s(wk​−wk−1​)c+dh0​,dhmax​,​ if wk​≤3 if wk​>3 if dhT,k​>dhmax​​(4) d h 0 dh_0 dh0​:初始高差阈值， s s s:坡度， c c c:网格尺寸， d h m a x dh_{max} dhmax​:最大高差阈值, k k k:迭代次数。 在城市地区，主要的非地面物体包括汽车、树木和建筑物。单个汽车和树木的尺寸远小于建筑物的尺寸，因此大多数汽车和树木通常在前几次迭代中被移除，而大型建筑物将在最后被移除。最大高差阈值可以设置为固定高度（例如，最低建筑高度），以确保识别建筑群。通常通过反复比较过滤和未过滤的数据来达到最佳效果。另一方面，山区的非地面物体主要是植被（树木）。不需要设置固定的最大高差阈值来移除树木，通常将其设置为研究区域内最大的高差。

1.计算 x , y x,y x,y最大值最小值 2.划分 m ∗ n m*n m∗n二维网格： m = m= m= floor [ ( max ⁡ ( y ) − min ⁡ ( y ) ) / c ] + 1 [(\max (y)-\min (y)) / c]+1 [(max(y)−min(y))/c]+1 and n = n= n= floor [ ( max ⁡ ( x ) − min ⁡ ( x ) ) / c ] + 1 [(\max (x)-\min (x)) / c]+1 [(max(x)−min(x))/c]+1 3.将点云坐标放进二维数组 A [ m , n ] A[m,n] A[m,n]（二维数组表示网格）中。遍历每个点，根据其x和y坐标，确定该点落在哪个网格中。如果同一网格单元中有多个点，选择高程最小的点。 4.使用最近邻法插值A中不包含任何点的网格。将这些插值网格的x和y坐标设置为零，以将它们与包含激光雷达点的网格单元区分开来。将A复制到B。用0初始化二维数组的元素。 5.用公式（1）或（2）计算 w k w_k wk​， w k < m a x − w i n d o w − s i z e w_kd h_T Z[j]−Zf​[j]>dhT​ then flag [ i , j ] = w k [i, j]=w_k [i,j]=wk​ 17. \quad\quad end for j j j loop 18. \quad end for i i i loop 19. \quad if ( d h T > d h max ⁡ ) \left(d h_T>d h_{\max }\right) (dhT​>dhmax​) 20. \quad d h T = d h max ⁡ \quad d h_T=d h_{\max } dhT​=dhmax​ 21. \quad else 22. \quad d h T = s ( w k − w k − 1 ) c + d h 0 \quad d h_T=s\left(w_k-w_{k-1}\right) c+d h_0 dhT​=s(wk​−wk−1​)c+dh0​ 23. end for window size loop 24. for i = 1 i=1 i=1 to m m m 25. \quad for j = 1 j=1 j=1 to n n n 26. \quad \quad if ( B [ i , j ] ( x ) > 0 (B[i, j](x)>0 (B[i,j](x)>0 and B [ i , j ] ( y ) > 0 ) B[i, j](y)>0) B[i,j](y)>0) 27. \quad \quad \quad if ( ( ( flag [ i , j ] = 0 ) [i, j]=0) [i,j]=0) 28. \quad \quad \quad \quad B [ i , j ] B[i, j] B[i,j] is a ground point 29. \quad \quad \quad else 30. \quad \quad \quad \quad B [ i , j ] B[i, j] B[i,j] is a nonground point 31. \quad end for j j j loop 32. end for i i i loop Erosion ⁡ ( Z , w k ) ‾ \underline{\operatorname{Erosion}\left(Z, w_k\right)} Erosion(Z,wk​)​ :

Dilation ( Z , w k ) \left(Z, w_k\right) (Z,wk​) :