CN109242786B - 一种适用于城市区域的自动化形态学滤波方法 - Google Patents
一种适用于城市区域的自动化形态学滤波方法 Download PDFInfo
- Publication number
- CN109242786B CN109242786B CN201810926973.6A CN201810926973A CN109242786B CN 109242786 B CN109242786 B CN 109242786B CN 201810926973 A CN201810926973 A CN 201810926973A CN 109242786 B CN109242786 B CN 109242786B
- Authority
- CN
- China
- Prior art keywords
- filtering
- morphological
- dsm
- window
- calculating
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000001914 filtration Methods 0.000 title claims abstract description 164
- 230000000877 morphologic effect Effects 0.000 title claims abstract description 66
- 238000000034 method Methods 0.000 title claims abstract description 51
- 230000008859 change Effects 0.000 claims abstract description 26
- 238000012886 linear function Methods 0.000 claims abstract description 6
- 230000003628 erosive effect Effects 0.000 claims description 7
- 238000004458 analytical method Methods 0.000 claims description 6
- 238000004891 communication Methods 0.000 claims description 5
- 230000010339 dilation Effects 0.000 claims description 4
- 230000008569 process Effects 0.000 claims description 4
- 230000009191 jumping Effects 0.000 claims description 3
- 238000005260 corrosion Methods 0.000 claims description 2
- 230000007797 corrosion Effects 0.000 claims description 2
- 230000004044 response Effects 0.000 claims description 2
- 238000012163 sequencing technique Methods 0.000 claims description 2
- ATJFFYVFTNAWJD-UHFFFAOYSA-N Tin Chemical compound [Sn] ATJFFYVFTNAWJD-UHFFFAOYSA-N 0.000 description 5
- 238000012952 Resampling Methods 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 4
- 238000001514 detection method Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 238000010835 comparative analysis Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 230000001788 irregular Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000001174 ascending effect Effects 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000004451 qualitative analysis Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/73—Deblurring; Sharpening
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10028—Range image; Depth image; 3D point clouds
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20024—Filtering details
- G06T2207/20032—Median filtering
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种适用于城市区域的自动化形态学滤波方法，包括：S1，将点云转换为深度图像，获取二维格网数据；S2，采用中值去噪法去除深度图像中的噪声数据；S3，设定滤波窗口尺度范围，采用形态学高帽运算对格网数据进行尺寸标记；S4，设定面积和粗糙度约束条件，对最大建筑物尺寸进行探测，同时确定最优滤波窗口；S5，依据形态学滤波结果计算各局部地形区域梯度变化，将滤波阈值设定为梯度变化的线性函数；S6，根据自适应变化的滤波阈值，按点基元逐点进行滤波实现。本发明能够实现最优滤波窗口的自动化确定，提升自动化程度。采用基于梯度变化的阈值设定，此阈值能够根据实际地形进行自适应更新，提高了对于复杂地形区域的鲁棒性及滤波精度。
Description
技术领域
本发明涉及地理空间信息系统技术领域，特别是涉及一种适用于城市区域的自动化形态学滤波方法。
背景技术
随着数字城市和智慧城市的快速发展，迫切需要我们对周围所处的城市地形环境有一个更准确的理解。近年来，机载LiDAR(Light Detection And Ranging)技术已经逐渐成为获取空间地理信息的一种遥感新技术。通过集成三维激光扫描仪、惯性导航系统以及全球定位系统，机载LiDAR能够获取来自地面目标物体的方位、距离以及表面特性。与传统的遥感方法相比，机载LiDAR受外界影响较小，能够24小时全天候地进行数据采集。此外，机载LiDAR系统发射的激光脉冲能够穿透植被，获取植被遮挡下的地形数据。此特点也有利于在城市区域进行道路信息采集。现如今，机载LiDAR技术已经广泛应用于城市规划的各个领域，例如三维城市建模、道路提取、电力线提取等。
为实现以上机载LiDAR后处理应用，一个非常重要的步骤便是要从点云数据中去除地物点而保留地形点，此过程通常称之为点云滤波。近年来，关于点云滤波的研究有很多。这些滤波算法按照滤波原理的不同，可以分为以下几类：基于坡度的点云滤波算法、基于曲面拟合的点云滤波算法、基于不规则三角网(TIN)的点云滤波算法以及基于形态学的滤波算法。
基于坡度的点云滤波算法通常都基于以下假设，即非地面点与地面点间的坡度值要大于地面点间的坡度值。因此，通过设置合理的坡度阈值可以将地面点与非地面点进行分离。但是此类算法的滤波结果过度依赖于坡度阈值的设定，如果阈值设置不合理就难以获得良好的滤波结果。基于曲面拟合的点云滤波算法是通过采用一定的插值拟合方法建立一地形曲面，通过设定一些滤波准则，例如点到曲面的距离，来渐进地将非地面点进行剔除。基于TIN的点云滤波算法通常采用渐进迭代的方式进行滤波实现。首先，获取地形各局部区域的最低点，然后利用这些种子点建立一个粗略的TIN。通过计算各个点到对应三角形的反复角以及反复距离，将符合阈值条件的点不断的加入到TIN 中。此过程一直迭代直到没有符合条件的点可加入TIN为止。
形态学滤波算法主要涉及到一系列的形态学运算，例如开运算、闭运算、膨胀运算、腐蚀运算。与以上三种滤波算法相比，形态学滤波算法原理简单、实现效率高，因此广泛应用于点云滤波中。对于形态学滤波算法，滤波窗口尺寸对滤波结果有直接的影响。小的滤波窗口只能滤除小的地物，而大的滤波窗口却容易削平地形。为解决此问题，不少形态学滤波算法采用渐进迭代地方式进行实现，滤波窗口由大到小变化，先滤波大的地物，然后在逐级地滤波小的地物。形态学滤波算法对于地形坡度变化较小的区域往往能够获得较好的滤波结果。
目前虽有一些适用于城市区域的形态学滤波方法，但还存在以下缺点：
(1)为获得较高的滤波精度往往需要不断地进行参数设置和阈值调节，尤其对于形态学滤波算法，最大滤波窗口需要大于测量区域建筑物的最大尺寸，但此尺寸是未知的，故只能通过反复试验探测最佳的最大滤波窗口以获得良好的滤波结果。复杂的参数设置不仅会降低算法的自动化程度而且也不利于经验缺乏的工作人员进行算法实现。
(2)难以平衡拒真误差和纳伪误差，即难以在有效滤除地物的同时保护地形细节，这使得在复杂的城市区域中难以获得良好的滤波精度。
发明内容
为此，本发的实施例提出一种适用于城市区域的自动化形态学滤波方法，解决现有技术中的自动化程度低、滤波精度低的问题。
一种适用于城市区域的自动化形态学滤波方法，包括以下步骤：
S1，将点云转换为深度图像，获取二维格网数据；
S2，采用中值去噪法去除深度图像中的噪声数据；
S3，设定滤波窗口尺度范围，采用形态学高帽运算对格网数据进行尺寸标记；
S4，设定面积和粗糙度约束条件，对最大建筑物尺寸进行探测，同时确定最优滤波窗口；
S5，依据形态学滤波结果计算各局部地形区域梯度变化，将滤波阈值设定为梯度变化的线性函数；
S6，根据自适应变化的滤波阈值，按点基元逐点进行滤波实现。
根据本发明提供的适用于城市区域的自动化形态学滤波方法，采用一系列的形态学高帽运算，对点云格网数据进行尺寸标记，通过设置面积和粗糙度的约束条件，能够自动地对测量区域内最大建筑物的尺寸进行探测，实现最优滤波窗口的自动化确定。为能有效地平衡拒真误差和纳伪误差，采用基于梯度变化的阈值设定，此阈值能够根据实际地形进行自适应更新，提高算法对于复杂地形区域的鲁棒性以及滤波精度。
因此，与现有技术相比，本发明提出的方法至少具有以下有益效果：
(1)提出了一种自动化形态学滤波方法，该方法不需要人为地进行参数设置和阈值调节，增强算法对各种复杂城市环境的适应性，解决了自动化程度低问题。
(2)通过采用一系列的形态学高帽运算以及相应的约束条件，实现对城市区域内最佳滤波窗口的自动化探测。
(3)依据形态学滤波结果，计算地形各局部区域的梯度变化，实现滤波阈值随地形变化的自适应更新，提升了滤波精度低。
另外，根据本发明上述的适用于城市区域的自动化形态学滤波方法，还可以具有如下附加的技术特征：
进一步地，所述步骤S2具体包括以下步骤：
S21，对二维格网数据DSM进行中值滤波，获取中值滤波结果DSMMedian；
S22，计算DSM和DSMMedian之间的高差；
S23，将高差大于5米的格网标记为噪声格网；
S24，将DSM中噪声格网的灰度值用DSMMedian中相应格网的灰度值进行更新。
进一步地，所述步骤S21中的中值滤波具体包括以下步骤：
S211，设计3×3的滤波模板，并用此模板对二维格网数据进行依次遍历；
S212，获取模板区域内各个格网的灰度值，并将这些灰度值进行顺序排序；
S213，取这些灰度值的中间数据作为模板中心格网新的灰度值。
进一步地，所述步骤S3具体包括以下步骤：
S31，初始化一二维格网数据Size_Flag，其大小与DSM相同；
S32，设定最大滤波窗口，并获取一系列按尺寸从小到大排序的滤波窗口集；
S33，依次采用这些滤波窗口用形态学高帽运算对DSM格网数据进行处理，计算不同滤波窗口下各个格网的高程变化响应值DSMtop(i,j)；
S34，计算不同滤波窗口对应的滤波阈值tre；
S35，在二维格网数据Size_Flag中，将DSMtop(i,j)大于阈值tre的格网标记为相对应滤波窗口的尺寸。
进一步地，所述步骤S33的形态学高帽运算具体包括以下步骤：
S331，对二维格网数据DSM进行进一步地行形态学腐蚀运算；
S332，进而对形态学腐蚀运算的结果进行形态学膨胀运算，获得开运算结果ODSM(i,j)；
S333，计算DSM(i,j)和ODSM(i,j)之间的差值，获得形态学高帽运算结果DSMtop(i,j)。
进一步地，所述步骤S4具体包括以下步骤：
S41，对二维格网数据Size_Flag进行四连通分析，获得连通结果 BW；
S42，设定面积约束阈值为T1；
S43，根据面积阈值计算出对应的粗糙度阈值；
S44，按从大到小的顺序，依次遍历BW中的各个连通部分，判断其是否符合面积和粗糙度的约束条件；
S45，如果不符合，则继续遍历；如果符合，则跳出循环遍历，此时连通部分所对应的滤波窗口即为最佳的滤波窗口。
进一步地，所述步骤S41中四连通分析具体包括以下步骤：
S411，依次遍历Size_Flag中各个格网的数值；
S412，如果在上下左右四个位置存在与该格网数值相同的格网，则将这两个格网标记为连通。
进一步地，所述步骤S5具体包括以下步骤：
S51，获取形态学滤波结果DTMMorph；
S52，计算获取DTMMorph各局部区域的梯度变化；
S53，将点基元滤波阈值Tr(k)表示为梯度变化的线性函数，实现滤波阈值随地形变化的自适应更新。
进一步地，所述步骤S6具体包括以下步骤：
计算各个点与DTMMorph之间的残差，将残差大于阈值Tr(k)的点标记为非地面点并进行剔除。
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实施例了解到。
附图说明
本发明实施例的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
图1是根据本发明提供的方法的流程图；
图2为最佳滤波窗口确定阶段伪代码流程图；
图3(a)为样本sample11的数字表面模型；
图3(b)为样本sample11的滤波后的数字地面模型；
图3(c)为样本sample11的准确的数字地面模型；
图4(a)为样本sample22的数字表面模型；
图4(b)为样本sample22的滤波后的数字地面模型；
图4(c)为样本sample22的准确的数字地面模型；
图5(a)为样本sample42的数字表面模型；
图5(b)为样本sample42的滤波后的数字地面模型；
图5(c)为样本sample42的准确的数字地面模型。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明一实施例提出的适用于城市区域的自动化形态学滤波方法，该方法能够自动对城市区域最大建筑物进行探测，从而实现最佳滤波窗口的确定，避免了人为的过多干预提高了算法的自动化程度。同时，该方法还能够依据地形变化，通过计算各局部地形区域的梯度变化，实现滤波阈值的自适应更新，提高了算法的鲁棒性以及对复杂城市环境的适应性。
本实施例中，涉及的方法的步骤流程图如图1所示，具体包括以下步骤：
S1，将点云转换为深度图像，获取二维格网数据；
S2，采用中值去噪法去除深度图像中的噪声数据；
S3，设定滤波窗口尺度范围，采用形态学高帽运算对格网数据进行尺寸标记；
S4，设定面积和粗糙度约束条件，对最大建筑物尺寸进行探测，同时确定最优滤波窗口；
S5，依据形态学滤波结果计算各局部地形区域梯度变化，将滤波阈值设定为梯度变化的线性函数；
S6，根据自适应变化的滤波阈值，按点基元逐点进行滤波实现。
其中，上述六个步骤构成三个阶段：Ⅰ深度图像去噪阶段(包括步骤S1和S2)；Ⅱ最佳滤波窗口确定阶段(包括步骤S3和S4)；Ⅲ滤波阈值自适应更新阶段(包括步骤S5和S6)。下面就这三个阶段进行具体说明。
Ⅰ深度图像去噪阶段
该阶段主要由以下两个环节组成，分别为将点云数据进行重采样生成二维格网数据DSM以及对该二维格网数据进行中值去噪。
为便于应用一些图像处理技巧，首先将点云数据进行重采样生成二维格网数据，以2.5维的深度图像形式进行表达。重采样的间隔一般设定为平均点间距，可采用下式进行估测：
式中，cellsize为格网尺寸，λ为平均点密度(个/平方米)。
格网尺寸确定后，就可对点云数据进行重采样，如果格网内有多于一个点存在，则该格网的灰度值为该格网内最低点的高程值；如果格网内没有无点云数据存在，则该格网的灰度值利用临近格网按照反比例加权法计算获取，公式表示如下。
式中，DSM(i,j)为需要插值计算的空白格网，DSM(i+u,j+v)为周围的临近格网的灰度值，di+u,j+v为点间距。
深度图像生成后，便可按中值滤波法对其进行去噪。中值滤波的模板选为3×3。通过计算中值滤波结果DSMMedian和原深度图像DSM之间的差值，将差值大于5米的格网标记为噪声格网，并用DSMMedian中相对应格网的灰度值对噪声格网进行灰度更新，获取去噪后的深度图像。
Ⅱ最佳滤波窗口确定阶段
最佳滤波窗口确定算法伪代码流程图如图2所示，具体包括以下步骤：
首先，设定滤波窗口的尺度范围Max_scales，生成一系列的按升序排列的滤波窗口集合，公式表示如下：
All_windows＝[1,2,...,Max_scales]
式中，Max_scales为一常数，通常设定为100。这是因为100×100的滤波窗口基本能实现对所有城市区域最大建筑物的探测。
然后，依次采用这些滤波窗口，对DSM进行形态学高帽运算，公式表示如下：
式中，DSM(i,j)为格网(i,j)处的灰度值，ODSM(i,j)为形态学开运算， EDSM(i,j)为形态学腐蚀运算，DDSM(i,j)为形态学膨胀运算。形态学开运算通常定义为先进行腐蚀运算，然后再对腐蚀运算结果进行膨胀运算。形态学腐蚀运算为取滤波窗口内灰度值的最小值作为该点新的灰度值，而形态学膨胀运算则是取滤波窗口内灰度值的最大值作为该点新的高程值，公式表示如下：
继而，依据滤波窗口尺寸大小的不同计算各个窗口所对应的阈值，公式表示如下：
式中，η和
分别设定为0.2和0.1。这是因为当滤波窗口取到最小值时，例如1×1，阈值计算为0.3，此阈值能够保证将地面点与低矮植被进行有效分离。
将形态学高帽运算结果DSMtop(i,j)与阈值tre进行比较，将大于阈值的格网标记为相对应的滤波窗口尺寸，获得尺寸标记结果Size_Flag。
将Size_Flag进行四连通分析，获得连通结果BW。为探测最大建筑物的尺寸，设定面积和粗糙度约束条件。这是因为建筑物的屋顶相较于树冠等地物来说，面积一般都比较大。此外，建筑物的屋顶一般都比较平坦，故粗糙度值相对较小。计算连通结果BW中各连通部分的粗糙度Roughbess和连通面积Area的方法，公式表示如下：
Area＝num*cellsize2
式中，num为BW中连通部分的总数，BWlabels为其中任一连通部分，cellsize格网尺寸。Mean表示计算均值，abs表示计算绝对值。
面积阈值T1设置为100，根据面积阈值可计算出相对应的粗糙度阈值T2，公式表示如下：
式中，
为计算连通部分中格网的数目，0.5为每个格网的容许误差。
依次遍历BW中的各个连通部分，如果存在满足面积和粗糙度约束条件的连通部分，则跳出循环并输出此时的滤波窗口为最佳滤波窗口。
Ⅲ滤波阈值自适应更新阶段
最佳滤波窗口确定后，便可进行形态学滤波运算，获取相对应的滤波结果。对于大多数的形态学滤波算法，点云数据都是重采样为二维格网。这样做虽然能够加快数据的处理效率，便于算法实现，但点云数据重采样往往会减低数据的精度，尤其是在一些断裂地形区域更容易产生误差。为解决此问题，可按点基元进行逐点滤波实现。各点的滤波阈值依据局部地形梯度变化自适应计算获取。公式表示如下：
Tr(k)＝Gradient(DTMMorph(i,j))+Δ
式中，Tr(k)为各点的滤波阈值，DTMMorph为形态学滤波结果，Δ为一常量，设置为0.5。Gradient表示计算局部地形区域的梯度，
和
分别表示计算地形局部区域x和y方向的梯度变化。各点的滤波阈值计算得到后，将各点到DTMMorph距离残差大于阈值Tr(k)的点标记为地物点并进行剔除。
为检验本发明所提出的滤波方法的有效性，选用ISPRS第三委员会发布的7组位于城市区域的专门用来检验滤波效果的点云数据进行实验。这7组样本数据由OptechALTM机载激光扫描仪获取，点间距为1-1.5m。这7组样本数据包含有多种复杂的城市地形特征，例如低矮植被、斜坡上建有房屋、不规则建筑物等。具体地形特征如表 1所示。这些滤波难点能够有效检验滤波算法的好坏。此外，这些样本数据都经过了人工编辑，地形点与地物点进行了有效的分离。因此，便于对滤波算法进行定量分析。
表1样本数据特征
采用四个精度指标对本发明的滤波效果进行评价，分别为Type Ⅰ误差、Type Ⅱ误差、总误差以及Kappa系数。Type Ⅰ误差又称拒真误差，指的是地面点误判为地物点的比例。Type Ⅱ误差又称纳伪误差，指的是地物点误判为地面点的比例。总误差指的是总的误判点的比例。Kappa系数指的是分类与完全随机的分类产生错误减少的比例。本发明针对这7组样本数据计算出来的这四类精度指标以及所得出的最佳滤波窗口的尺寸如表2所示。
表2精度计算结果以及最佳滤波窗口
为了更加直观地验证本发明滤波算法的有效性，分别选用 sample11、sample22和sample42进行定性分析。这三组样本的滤波前的数字表面模型(a)、滤波后的数字地面模型(b)以及准确的数字地面模型(c)分别如图3、4、5所示。从图中可以看出，本发明的滤波结果(图(b))与准确的滤波结果(图(c))十分接近。结合表2，本发明提出的滤波方法能够有效地对各种复杂的城市环境实现点云滤波，三类误差的平均值均小于5％，平均Kappa系数高于90％。而且本发明的滤波方法可以取得相对均衡的Type Ⅰ误差(平均值为 4.70％)和TypeⅡ误差(平均值为3.62％)，由此可以看出本发明在有效去除地物的同时能够尽可能地保护地形细节。
为进一步说明本发明的滤波效果，选用现有技术中的其他7种滤波方法(对照方法1～7)对总误差和Kappa系数进行对比分析，分析结果如表3和4所示。从对比结果可以看出，本发明方法能够取得最小的平均总误差以及最大的平均Kappa系数。由此可以说明，本发明方法能够取得更好的滤波效果。此外，本研究还计算了七组样本数据滤波结果的最小值(min)、最大值(max)以及标准差(std)。从这三组统计数据可以看出，本发明方法能够取得最小的标准差(总误差的标准差为2.39，Kappa系数的标准差为4.98)，由此可进一步说明本发明方法对于各种样本数据都能取得较好的滤波结果，因此本发明方法鲁棒性更好，能够适应多种复杂的城市环境。
表3 7组样本数据总误差对比分析
表4 7组样本数据Kappa对比分析(后两种方法未提供Kappa信息)
与现有技术相比，本实施例提出的方法至少具有以下有益效果：
(1)提出了一种自动化形态学滤波方法，该方法不需要人为地进行参数设置和阈值调节，增强算法对各种复杂城市环境的适应性，解决了自动化程度低问题。
(2)通过采用一系列的形态学高帽运算以及相应的约束条件，实现对城市区域内最佳滤波窗口的自动化探测。
(3)依据形态学滤波结果，计算地形各局部区域的梯度变化，实现滤波阈值随地形变化的自适应更新，提升了滤波精度低。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。
Claims (9)
1.一种适用于城市区域的自动化形态学滤波方法，其特征在于，包括以下步骤：
S1，将点云转换为深度图像，获取二维格网数据；
S2，采用中值去噪法去除深度图像中的噪声数据；
S3，设定滤波窗口尺度范围，采用形态学高帽运算对格网数据进行尺寸标记；
S4，设定面积和粗糙度约束条件，对最大建筑物尺寸进行探测，同时确定最佳 滤波窗口；
S5，通过最佳滤波窗口进行形态学滤波运算，获取形态学滤波结果，依据形态学滤波结果计算各局部地形区域梯度变化，将滤波阈值设定为梯度变化的线性函数；
S6，根据自适应变化的滤波阈值，按点基元逐点进行滤波实现。
2.根据权利要求1所述的适用于城市区域的自动化形态学滤波方法，其特征在于，所述步骤S2具体包括以下步骤：
S21，对二维格网数据DSM进行中值滤波，获取中值滤波结果DSMMedian；
S22，计算DSM和DSMMedian之间的高差；
S23，将高差大于5米的格网标记为噪声格网；
S24，将DSM中噪声格网的灰度值用DSMMedian中相应格网的灰度值进行更新。
3.根据权利要求2所述的适用于城市区域的自动化形态学滤波方法，其特征在于，所述步骤S21中的中值滤波具体包括以下步骤：
S211，设计3×3的滤波模板，并用此模板对二维格网数据进行依次遍历；
S212，获取模板区域内各个格网的灰度值，并将这些灰度值进行顺序排序；
S213，取这些灰度值的中间数据作为模板中心格网新的灰度值。
4.根据权利要求3所述的适用于城市区域的自动化形态学滤波方法，其特征在于，所述步骤S3具体包括以下步骤：
S31，初始化一二维格网数据Size_Flag，其大小与DSM相同；
S32，设定最大滤波窗口，并获取一系列按尺寸从小到大排序的滤波窗口集；
S33，依次采用这些滤波窗口用形态学高帽运算对DSM格网数据进行处理，计算不同滤波窗口下各个格网的高程变化响应值DSMtop(i,j)；
S34，计算不同滤波窗口对应的滤波阈值tre；
S35，在二维格网数据Size_Flag中，将DSMtop(i,j)大于阈值tre的格网标记为相对应滤波窗口的尺寸。
5.根据权利要求4所述的适用于城市区域的自动化形态学滤波方法，其特征在于，所述步骤S33的形态学高帽运算具体包括以下步骤：
S331，对二维格网数据DSM进行进一步地行形态学腐蚀运算；
S332，进而对形态学腐蚀运算的结果进行形态学膨胀运算，获得开运算结果ODSM(i,j)；
S333，计算DSM(i,j)和ODSM(i,j)之间的差值，获得形态学高帽运算结果DSMtop(i,j)。
6.根据权利要求5所述的适用于城市区域的自动化形态学滤波方法，其特征在于，所述步骤S4具体包括以下步骤：
S41，对二维格网数据Size_Flag进行四连通分析，获得连通结果BW；
S42，设定面积约束阈值为T1；
S43，根据面积阈值计算出对应的粗糙度阈值；
S44，按从大到小的顺序，依次遍历BW中的各个连通部分，判断其是否符合面积和粗糙度的约束条件；
S45，如果不符合，则继续遍历；如果符合，则跳出循环遍历，此时连通部分所对应的滤波窗口即为最佳的滤波窗口。
7.根据权利要求6所述的适用于城市区域的自动化形态学滤波方法，其特征在于，所述步骤S41中四连通分析具体包括以下步骤：
S411，依次遍历Size_Flag中各个格网的数值；
S412，如果在上下左右四个位置存在与该格网数值相同的格网，则将这两个格网标记为连通。
8.根据权利要求7所述的适用于城市区域的自动化形态学滤波方法，其特征在于，所述步骤S5具体包括以下步骤：
S51，通过最佳滤波窗口进行形态学滤波运算，获取形态学滤波结果DTMMorph；
S52，计算获取DTMMorph各局部区域的梯度变化；
S53，将点基元滤波阈值Tr(k)表示为梯度变化的线性函数，实现滤波阈值随地形变化的自适应更新。
9.根据权利要求8所述的适用于城市区域的自动化形态学滤波方法，其特征在于，所述步骤S6具体包括以下步骤：
计算各个点与DTMMorph之间的残差，将残差大于阈值Tr(k)的点标记为非地面点并进行剔除。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810926973.6A CN109242786B (zh) | 2018-08-15 | 2018-08-15 | 一种适用于城市区域的自动化形态学滤波方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810926973.6A CN109242786B (zh) | 2018-08-15 | 2018-08-15 | 一种适用于城市区域的自动化形态学滤波方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109242786A CN109242786A (zh) | 2019-01-18 |
| CN109242786B true CN109242786B (zh) | 2020-09-22 |
Family
ID=65071355
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810926973.6A Active CN109242786B (zh) | 2018-08-15 | 2018-08-15 | 一种适用于城市区域的自动化形态学滤波方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109242786B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110119438B (zh) * | 2019-04-23 | 2021-05-07 | 东华理工大学 | 基于主动学习的机载LiDAR点云滤波方法 |
| CN112734819B (zh) * | 2021-01-14 | 2022-04-19 | 长光卫星技术股份有限公司 | 一种适用于高分辨率遥感卫星dsm的地表滤波方法 |
| CN114332631B (zh) * | 2022-01-12 | 2023-04-18 | 中铁二院工程集团有限责任公司 | 一种适用于山区危岩落石的LiDAR点云数据提取方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102520401A (zh) * | 2011-12-21 | 2012-06-27 | 南京大学 | 一种基于LiDAR数据的建筑物区域提取方法 |
| CN103745436A (zh) * | 2013-12-23 | 2014-04-23 | 西安电子科技大学 | 基于区域预测的LiDAR点云数据形态学滤波方法 |
| CN105787921A (zh) * | 2015-08-19 | 2016-07-20 | 南京大学 | 利用机载LiDAR数据重建大型复杂立交桥三维模型的方法 |
| WO2017120897A1 (zh) * | 2016-01-15 | 2017-07-20 | 武汉武大卓越科技有限责任公司 | 基于线扫描三维点云的物体表面变形特征提取方法 |
-
2018
- 2018-08-15 CN CN201810926973.6A patent/CN109242786B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102520401A (zh) * | 2011-12-21 | 2012-06-27 | 南京大学 | 一种基于LiDAR数据的建筑物区域提取方法 |
| CN103745436A (zh) * | 2013-12-23 | 2014-04-23 | 西安电子科技大学 | 基于区域预测的LiDAR点云数据形态学滤波方法 |
| CN105787921A (zh) * | 2015-08-19 | 2016-07-20 | 南京大学 | 利用机载LiDAR数据重建大型复杂立交桥三维模型的方法 |
| WO2017120897A1 (zh) * | 2016-01-15 | 2017-07-20 | 武汉武大卓越科技有限责任公司 | 基于线扫描三维点云的物体表面变形特征提取方法 |
Non-Patent Citations (3)
| Title |
|---|
| A progressive morphological filter for point cloud extracted from UAV images;Qiuling Wang等;《2014 IEEE Geoscience and Remote Sensing Symposium》;20140718;第2023-2026页 * |
| 决策树约束的建筑点云提取方法;雷钊 等;《激光与光电子学进展》;20180627;第1-7页 * |
| 基于改进窗口尺寸的LiDAR点云数据滤波;孙涛 等;《江西科学》;20180210;第36卷(第1期);第150-155页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109242786A (zh) | 2019-01-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Tóvári et al. | Segmentation based robust interpolation-a new approach to laser data filtering | |
| CN105844629B (zh) | 一种大场景城市建筑物立面点云自动分割方法 | |
| JP4545219B1 (ja) | 地形画像を用いた地形変化の解析方法及びそのプログラム | |
| CN109242786B (zh) | 一种适用于城市区域的自动化形态学滤波方法 | |
| CN111598780B (zh) | 一种适用于机载LiDAR点云的地形自适应插值滤波方法 | |
| CN111738945B (zh) | 一种基于矿山的点云数据预处理方法 | |
| CN111007531A (zh) | 一种基于激光点云数据的道路边沿检测方法 | |
| CN112945196B (zh) | 一种基于点云数据的露天矿台阶线提取和边坡监测的方法 | |
| CN112907744B (zh) | 数字高程模型的构建方法、装置、设备及存储介质 | |
| Chen et al. | A mathematical morphology-based multi-level filter of LiDAR data for generating DTMs | |
| CN115512247A (zh) | 基于图像多参数提取的区域建筑损伤等级评定方法 | |
| CN112734819B (zh) | 一种适用于高分辨率遥感卫星dsm的地表滤波方法 | |
| Brzank et al. | Automated extraction of pair wise structure lines using airborne laserscanner data in coastal areas | |
| CN111080536A (zh) | 一种机载激光雷达点云的自适应滤波方法 | |
| CN111583406A (zh) | 杆塔脚基点坐标计算方法、装置及终端设备 | |
| CN115294302A (zh) | 一种基于断裂线约束的机载点云快速滤波方法 | |
| CN112381029B (zh) | 一种基于欧氏距离的机载LiDAR数据建筑物提取方法 | |
| Chang et al. | Bare-earth extraction from airborne LiDAR data based on segmentation modeling and iterative surface corrections | |
| Bucksch et al. | Skeleton-based botanic tree diameter estimation from dense LiDAR data | |
| CN114972358B (zh) | 一种基于人工智能的城市测绘激光点云偏移检测方法 | |
| Nizar et al. | Reconstruction of buildings from airborne laser scanning data | |
| CN118072029B (zh) | 一种改进泰森多边形约束的车载点云单木分割方法及系统 | |
| Štroner et al. | Multi-View and Shift Rasterization Algorithm (MVSR) for Effective Identification of Ground in Dense Point Clouds | |
| Lu et al. | LiDAR Data Filtering Based on the Improved Window Size of Progressive Mathematical Morphology | |
| Ojiako et al. | REMOVAL OF NOISE FROM SHUTTLE RADAR TOPOGRAPHIC MISSION OVER VEGETATION AND BUILT UP AREA FOR LARGE SCALE TOPGRAPHIC MAPPING |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |