CN112308966A - 基于多级曲率约束的高斯混合模型点云滤波方法 - Google Patents
基于多级曲率约束的高斯混合模型点云滤波方法 Download PDFInfo
- Publication number
- CN112308966A CN112308966A CN202011253108.3A CN202011253108A CN112308966A CN 112308966 A CN112308966 A CN 112308966A CN 202011253108 A CN202011253108 A CN 202011253108A CN 112308966 A CN112308966 A CN 112308966A
- Authority
- CN
- China
- Prior art keywords
- point
- ground
- points
- point cloud
- classified
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 69
- 238000001914 filtration Methods 0.000 title claims abstract description 46
- 239000000203 mixture Substances 0.000 title claims abstract description 37
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 31
- 238000004364 calculation method Methods 0.000 claims abstract description 18
- 230000008569 process Effects 0.000 claims abstract description 14
- 238000010276 construction Methods 0.000 claims description 8
- 239000011159 matrix material Substances 0.000 claims description 4
- 238000011160 research Methods 0.000 claims description 4
- 239000000126 substance Substances 0.000 claims description 4
- 238000005457 optimization Methods 0.000 claims description 3
- 230000002159 abnormal effect Effects 0.000 abstract description 2
- 238000013459 approach Methods 0.000 abstract description 2
- 238000005259 measurement Methods 0.000 description 9
- 238000010586 diagram Methods 0.000 description 8
- 238000012545 processing Methods 0.000 description 8
- 238000012360 testing method Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 3
- 230000011218 segmentation Effects 0.000 description 3
- 238000013507 mapping Methods 0.000 description 2
- 230000000877 morphologic effect Effects 0.000 description 2
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011439 discrete element method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000007246 mechanism Effects 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
- 238000003908 quality control method Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
- G01S17/02—Systems using the reflection of electromagnetic waves other than radio waves
- G01S17/06—Systems determining position data of a target
- G01S17/42—Simultaneous measurement of distance and other co-ordinates
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
- G01S17/88—Lidar systems specially adapted for specific applications
- G01S17/89—Lidar systems specially adapted for specific applications for mapping or imaging
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
- G01S7/4802—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
- G06T17/205—Re-meshing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Theoretical Computer Science (AREA)
- Software Systems (AREA)
- Geometry (AREA)
- Computer Networks & Wireless Communication (AREA)
- Electromagnetism (AREA)
- Radar, Positioning & Navigation (AREA)
- Computer Graphics (AREA)
- Processing Or Creating Images (AREA)
- Image Processing (AREA)
Abstract
本发明提出了一种基于多级曲率约束的高斯混合模型点云滤波方法，该方法主要包括三个步骤：(1)根据大数定律，首先计算每个测量点到其K近邻点的平均距离，然后基于正常测量值不会大于M倍标准差的原则，剔除点云中存在的各种异常测量值。(2)其次，对由局部最小值获得的地面种子点，利用薄板样条函数(TPS)进行插值，然后剔除与插值曲面高差大于曲率阈值的点。该过程是通过不断减少的插值格网单元与逐渐增加曲率阈值来迭代逼近一个参考曲面。(3)最后，计算每一个测量值的相对高差，应用二元高斯混合模型(GMM)进行地物与地面的分类，再利用期望最大化算法(EM)得到模型参数的最优解，并根据计算得到后验概率判断该点是否属于地物完成点云滤波。
Description
技术领域
本发明属于激光雷达应用、计算机视觉和地理信息系统领域，主要包括机载LiDAR点云的地面滤波过程中的点云去噪、曲面格网构建以及根据后验概率确定点的类别等方面。
背景技术
机载LiDAR系统通过测量激光回波的时间来得到目标至仪器中心的距离。只要给定地理参考系，系统能够快速得到高精度的点云数据。与传统由摄影测量得到点云的方法比较， LiDAR系统的优势在于较短的处理时间以及更高的测量精度。当前，LiDAR技术已经被推广到多个领域，如数字地形模型构建、三维建筑模型构建、路面提取以及森林站参数的估计等等。在众多的应用中，将地面测量与非地面测量在原始点云中区分开(滤波)是一项必不可少的工作。同样，精确的数字高程模型(DEM)构建也需要预先将非地面点移除，进而由地学插值方法计算得到栅格格网。
机载LiDAR技术革新了获取高精度DEM的方式，这主要是因为其获取点云效率高、精度高以及价格低廉。相较于数字摄影测量技术，机载LiDAR提供了一种对大范围区域的高精度地形测绘方案，许多研究者将注意力集中于如何使用LiDAR点云数据获取地形信息。为了获取DEM，必须要进行的步骤是过滤掉非地面点云，包括来自桥梁、建筑、植被等的测量值。值得注意的是，虽然机载LiDAR技术能够在较短的时间内获得庞大的数据量，但是其滤波以及手动质量控制却要占到整个后处理时间的60％-80％。因此，如何减少点云数据处理算法的时间复杂度并减少其参数调整时间是当前研究工作的主要问题。
在过去的二十年间，有数十种地面滤波算法被提出来。然而，从处理时间以及处理精度来看，滤波依旧是点云数据处理的一大挑战。Sithole和Vosselman在2004年对8种经典的滤波算法于农村以及城市试验区数据开展了综合性的比较实验工作。他们指出所有的算法在较为平坦的地形上都能够表现出色。然而，当地形变得陡峭并伴随有复杂的地物(如桥梁、低矮灌木等)，滤波器将会失效而表现得不尽人意。总体而言，基于曲面的滤波器在大多数场景的表现大多数较好。这主要是因为该类型滤波器相较于其他类型滤波器利用了更多的点云数据上下文信息。自此，他们选用的位于农村及城市的15块测试样本被用作测试地面滤波算法的标准数据集。
机载LiDAR点云的地面滤波算法主要可以分为四大类：基于坡度的滤波器，基于数学形态学的滤波器，基于曲面插值的滤波器，基于分割策略的滤波器。其中基于数学形态学以及基于曲面插值的滤波器是广泛使用的方法。特别地，商用软件TerraScan集成了经典渐进三角加密滤波(PTD)算法，能够以相对高的精度得到地形测量。但是，该算法的一大不足就是需要复杂的参数调整，且因三角面片会在特殊地形上变形严重而无法满足任意地形的地面滤波精度要求。而基于坡度的滤波器多假设地形点测量间的坡度将小于非地形点，通过给定坡度阈值滤除非地形点，由于假设条件大多适用于较为平坦的地形，而大多数的地形情况是复杂多变的，因此限制了该类方法在实际工作中的应用。基于分割策略的滤波器的基本思想是先聚类后分类，通过同质性原则采用聚类算法如区域增长方法进行测量点的合并，然后给定角度或者距离阈值对各个分割段分类为地面段与非地面段，从而达到滤波目的，但该方法通常需要选择合适的聚类方法保证分割段内部均质性，且伴随着复杂参数设置问题，在滤波工作中较少使用。
现今的大多数滤波算法都需要根据具体的试验区域进行复杂的参数调整才能达到要求的精度。然而，这对于许多测绘业界工作者而言这并不是十分友好，因为他们需要对算法机制有着一定的了解且需要大量的实验调整参数。举例而言，当使用基于数学形态学的滤波器时，算法往往假定坡度s为一恒定的常数。但这个假设在实际工作中无法得到满足，因为受到地表内、外力作用下，地形坡度总是不断变化的。使用恒定的s值可能会导致，高度阈值dT计算错误，最终产生大量的遗漏误差或者委托误差，这在后续DEM构建中将导致灾难性错误。
发明内容
本发明的目的在于针对现有地面滤波算法存在的复杂参数调整、滤波结果需要后期交互验证修正的不足，提出一种新型的基于多级曲率约束的高斯混合模型点云滤波方法。
基于多级曲率约束的高斯混合模型点云滤波方法，包括以下步骤：
步骤一，输入机载LiDAR点云，根据大数定律，首先计算每个点云到其K近邻点的平均距离，然后基于大于M倍标准差为粗差的原则，剔除不满足条件的点云，得到待处理点云；
步骤二，对待处理点云所在区域建立网格索引，并由局部最小值方法获得的地面种子点，利用薄板样条函数(TPS)进行插值，剔除与插值曲面高差大于曲率阈值的点，得到待分类地面点；
步骤三，计算每一个待分类地面点的相对高差ΔH，应用二元高斯混合模型(GMM)进行地物与地面的分类，使用期望最大化算法(EM)得到模型参数的最优解，并根据计算得到后验概率判断该点是否属于地物完成点云滤波。
所述步骤一的具体过程为：
(1)输入机载LiDAR点云，并根据研究区域点云分布，由用户给定离群点去噪算法(SOR) 参数，SOR参数包括近邻点个数K以及标准差乘系数M；
(2)根据大数定律，遍历每个点云计算到其K近邻点的平均距离，再基于大于M倍标准差为粗差的原则，剔除不满足
的点云，得到待处理点云；其中，
为第i个点云到其K近邻点的平均距离，KNN平均距离
标准差
n为输入的机载LiDAR点云数。
所述步骤二的具体过程为：
(1)先对步骤一得到的待处理点云，由局部最小值法获取初始地面种子点集Gseed；给定格网大小gw，构建网格索引，根据网格行列号(m,n)计算出网格单元中心的真实平面坐标值(x, y)，再由最小化薄板样条能量函数Ef的优化原则，选择待插值格网点的N个近邻地面种子点作为控制点；
(2)令控制点权重满足插值函数f(x,y)存在二阶偏导数约束条件，根据拉格朗日乘数法，解算出插值函数f(x,y)中的未知系数，进而插值出网格单元中心高程值；
(3)给定曲率阈值参数t，由gw和t构建曲面格网金字塔R1，并使用一个3×3的均值移动窗口经过每个待处理点云所在格网，得到一个平滑后的曲面格网Rnew，根据Rnew高程值znew(m,n)进而计算平均曲率c(m,n)＝znew(m,n)+t；遍历待处理点云，若第j个待处理点云Pj(xj,yj,zj)的高程值zj大于等于c(m,n)，则将其加入非地面点，否则将其作为第一待分类地面点进入(4)；
(4)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔R2，并根据(3)中的方法遍历第一待分类地面点筛选出第二待分类地面点，进入(5)；
(5)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔R3，并根据(3)中的方法遍历第二待分类地面点筛选出第三待分类地面点，第一至第三待分类地面点统称待分类地面点。
其中，(1)中根据下式换算真实平面坐标值(x，y)：
式中，(xmin，ymin)为输入的机载LiDAR点云的最小边界盒的x、y方向上最小坐标值。
其中，(2)中插值函数：
式中，a1、a2、a3为反映插值曲面总体趋势的系数；ws表示与第s个控制点cs相关的权重系数；
为与薄板样条函数(TPS)相关的径向基核函数，(xs，ys)为cs的x、y方向坐标值，N为控制点个数；
f(x,y)的二阶偏导数约束为：
式中，zs为cs在z方向的坐标值；
根据拉格朗日乘数法，利用以下线性方程组求解f(x,y)中的未知系数：
式中，矩阵U的第I行第J列元素
(xI，yI)和(xJ， yJ)分别为第I控制点cI和第J个控制点cJ的x、y方向坐标值；
最后，待插值单元高程值z的计算公式为：z＝f(x,y)。
所述步骤三的具体过程为：
(1)遍历所有的待分类点，计算每一个待分类点的相对高程，将ΔH代入二元高斯混合模型(GMM)中，计算高斯混合模型对每个待分类点的相对高程的响应度，应用期望最大化算法(EM)迭代更新计算各分模型参数值αl,μl,σl 2,l＝1,2；其中，第k个待分类点qk的相对高程ΔHk＝H–Hgrid，H为待分类点的高程值，Hgrid为待分类点所在R3的曲面格网单元的高程，αl为模型的混合系数且大于等于0，μl是qk的相对高程的平均值；σl是qk的相对高程的标准差；
(2)由于假设地面点、非地面点的相对高程均符合高斯分布，所以根据计算出的各分模型高斯参数值，得到包含隐变量的条件概率估计值，进而使用贝叶斯公式及全概率公式求解用于区分地面点的后验概率P(G|qk)，若后验概率P(G|qk)大于等于0.5，则认为第k个待分类点qk属于地面点，并以之参与数字高程模型DEM构建；若后验概率P(G|qk)小于0.5，则认为qk属于非地面点，并舍弃该点；
其中，根据下式更新计算各分模型参数值αl,μl,σl 2,l＝1,2：
式中，ΔHk为待分类地面点qk的相对高程；P(G|qk)为待分类点qk属于地面点的后验概率； P(NG|qk)为待分类点qk属于非地面点的后验概率；Nc是待处理点云总数量。
其中，(2)中根据下式计算待分类点qk属于地面点的后验概率P(G|qk)及属于非地面点的后验概率P(NG|qk)：
式中，P(G)为待分类点qk属于地面点的先验概率；P(NG)为待分类点qk属于非地面点的先验概率；P(qk)由全概率公式计算，P(qk)＝P(qk|G)P(G)+P(qk|NG)P(NG)。
每个待分类点qk的条件概率：
式中，αl为模型的混合系数且大于等于0；ΔHk是待处理点qk的相对高程；φ()是高斯密度函数，θl＝(μl,σl)是其参数，μl是待分类点qk相对高程的平均值；σl是待分类点qk相对高程的标准差。
应用上式计算地面点的条件概率P(qk|G)，仅用到一组参数(α1,μ1,σ1)：
P(qk|G)＝α1φ(ΔHk|μ1,σ1)
同样地，计算非地面点的条件概率P(qk|NG)时用另一组参数(α2,μ2,σ2)：
P(qk|NG)＝α2φ(ΔHk|μ2,σ2)
其中高斯密度函数为：
本发明的基于多级曲率约束的高斯混合模型点云滤波方法是一种自适应的滤波方法，其结合了统计方法与曲面插值方法的优点，算法给出了较为通用的一组参数，能够方便有滤波工作需要的用户使用。本发明的滤波方法采用了曲面插值策略，通过薄板样条函数(TPS)构建曲面金字塔，避免了因局部地形坡度变化而导致的地形点误分类，通过曲率约束保证了地物点能够被有效区分并被剔除；此外，对每一个测量点进行相对高差的计算，能够有效地克服陡峭地形如悬崖上测量点绝对高程高于平坦地形上的树木测量点绝对高程而导致的误分类问题，进而使用二元高斯混合模型(GMM)通过概率判别的方式完成滤波，模型的求解由期望最大化算法(EM)自动地寻优、避免了过多主观判断、及减少了大量的人工调参时间。
附图说明
图1是本发明实施例的方法框架图；
图2是本发明实施例的统计局外点移除示意图，其中，(a)是移除前示意图，(b)是确定局外点集，(c)是移除后示意图；
图3是本发明实施例的多级曲率约束的曲面格网示意图，其中，(a)是选用ISPRS第三委员会提供的测试样本点云数据所构建的数字表面模型(DSM)，(b)是本发明实施例构建的曲面格网金字塔中第一级曲面格网R1，(c)是本发明实施例构建的曲面格网金字塔中第二级曲面格网R2，(d)是本发明实施例构建的曲面格网金字塔中第三级曲面格网R3；
图4是本发明实施例的一个实测点云滤波结果图，其中，(a)是使用原始LiDAR点云构建的数字表面模型(DSM)，(b)是使用本发明的基于多级曲率约束的高斯混合模型点云滤波方法得到的数字高程模型(DEM)；(c)是使用经典数学形态学滤波方法得到的数字高程模型(DEM)。
具体实施方式
以下结合具体实施例，并参照附图，对本发明作进一步详细说明。
本发明根据大数定律，首先计算每个测量点到其K近邻点的平均距离，然后基于正常测量值不会大于M倍标准差的原则，剔除点云中存在的各种异常测量值；其次，对由局部最小值获得的地面种子点，利用薄板样条函数(TPS)进行插值，然后剔除与插值曲面高差大于曲率阈值的点。该过程是通过不断减少的插值格网单元与逐渐增加曲率阈值来迭代逼近一个参考曲面；最后，计算每一个测量值的相对高差，应用二元高斯混合模型(GMM)进行地物与地面的分类，再利用期望最大化算法(EM)得到模型参数的最优解，并根据计算得到后验概率判断该点是否属于地物完成点云滤波，以有效地减少人工调参时间，提升机载LiDAR点云滤波精度与效率。
如附图1所示，本发明的基于多级曲率约束的高斯混合模型点云滤波方法包括三部分：①根据大数定律确定点云中的待处理点；②基于薄板样条插值的多级曲率约束；③相对高程计算与高斯混合模型分类。具体实施步骤为：
第一步：根据大数定律确定点云中的待处理点。
机载LiDAR点云数据中存在两类噪声点：高位噪声点(high outliers)、低位噪声点(low outliers)，利用本发明中的KNN邻域搜索可以检测并剔除这两类噪声。
根据大数定律确定点云中的待处理点的示意图如附图2中的(a)至(c)所示，其具体流程如下：
(1)输入机载LiDAR点云，并根据研究区域点云分布，由用户给定离群点去噪算法(SOR) 参数，SOR参数包括近邻点个数K以及标准差乘系数M；
(2)根据大数定律，遍历每个点云计算到其K近邻点的平均距离，再基于大于M倍标准差为粗差的原则，剔除不满足
的点云，得到待处理点云；其中，
为第i个点云到其K近邻点的平均距离，KNN平均距离
标准差
n为输入的机载LiDAR点云数。
假设待计算点云qi(xi，yi，zi)到其最近的K个邻居点pm(xm，ym，zm)，(m＝1,2,...,K)的距离为dim，则其计算公式如下：
则其KNN平均距离
的计算可以由下式给出：
式中，K的值一般由用户指定。
再根据
可以计算KNN平均距离dmean及标准差stddev，可以由统计学计算公式给出：
式中，n为输入机载LiDAR点云的总数量。
因此，令标准差乘系数M＝3，则待处理点云qi的筛选判断条件如下：
若点满足上式条件约束，则将其加入待处理点集Q中；若点不满足上式条件，则将其加入局外点集O并舍弃之。
第二步：基于薄板样条插值的多级曲率约束。
基于薄板样条插值的多级曲率约束示意图如附图3中的(a)至(d)所示，其过程为：
(1)先对步骤一得到的待处理点云，由局部最小值法获取初始地面种子点集Gseed，给定格网大小gw，构建网格索引，根据网格行列号(m,n)计算出网格单元中心的真实平面坐标值(x, y)，再由最小化薄板样条能量函数Ef的优化原则，选择待插值格网点的N个近邻地面种子点作为控制点；
假设待插值格网的行列号为(m,n)，根据下式换算成相应的真实平面坐标值(x，y)：
式中，gw为格网大小；(xmin，ymin)为原始点云最小边界盒的x,y方向上最小坐标值。
给定点三维空间中的N＝12个控制点cs(xs，ys，zs)，薄板样条函数(TPS)通过构造函数f(x, y)穿过所有的控制点cs并最小化弯曲能量函数Ef来完成插值。在二维空间中，Ef被定义为二阶偏导数平方的积分，由下式给出：
式中，f为插值函数f(x,y)；(dx,dy)为f(x,y)在x,y方向的微分。
(2)令控制点权重满足插值函数f(x,y)存在二阶偏导数约束条件，根据拉格朗日乘数法，解算出插值函数f(x,y)中的未知系数，进而插值出网格单元中心高程值；
以待插值单元gidx为例，由公式(6)导出其真实平面坐标值(x，y)，其高程值z按下面的方式计算：
首先，需要计算f(x,y)中的未知参数，其计算公式由下式给出：
式中，a1,a2,a3为反映插值曲面总体趋势的系数；ws表示与控制点cs相关的权重系数；
为与薄板样条函数(TPS)相关的径向基核函数，(xs，ys)为控制点cs的x，y方向坐标值；同时，为了满足f(x,y)存在二阶偏导数约束，需要满足下式条件。
然后，根据拉格朗日乘数法，利用线性方程组求解f(x,y)中的未知系数，其计算公式如下：
式中，U为由
定义的矩阵，N为控制点个数， (xI，yI),(xJ，yJ)为控制点的x，y方向坐标值；而w为由ws构成的列向量；a为由a1,a2,a3构成的列向量；相应地，P和v可以由下式表出：
其中，矩阵中的元素来源于一系列控制点集cs(xs，ys，zs)，(s＝1,2,…,N)，N为控制点个数。
最后，待插值单元gidx高程值z的计算公式为：
z＝f(x,y) (12)
(3)给定曲率阈值参数t，由gw和t构建曲面格网金字塔R1，并使用一个3×3的均值移动窗口经过每个待处理点云所在格网，得到一个平滑后的曲面格网Rnew，根据Rnew高程值znew(m,n)进而计算平均曲率c(m,n)＝znew(m,n)+t；遍历待处理点云，若点Pj(xj,yj,zj)的高程值zj大于等于c(m,n)，则将其加入非地面点，否则将其作为第一待分类地面点进入(4)。
由局部最小值法获取初始地面种子点集Gseed，然后给定格网大小gw与曲率阈值t，由公式z＝f(x,y)计算出插值曲面格网R1，然后遍历待处理点qj，并由公式(6)计算其所在的格网单元行列号(m,n)，采用图像卷积操作，利用一个3×3的均值移动窗口处理该格网单元，计算得到一个平滑曲面Rnew，根据Rnew高程值znew(m,n)，并按下式计算平均曲率c(m,n)：
c(m,n)＝znew(m,n)+t (13)
式中，znew(m,n)表示行列号(m,n)所在的平滑曲面Rnew格网单元；t为曲率阈值参数。
判断待处理点qj(xj,yj,zj)的高程值zj是否满足条件zj<c(m,n)，若是则将其加入待分类地面点集G；否则，将其加入非地面点集NG。
(4)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔R2，并根据(3)中的方法遍历第一待分类地面点筛选出第二待分类地面点，进入(5)；
(5)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔R3，并根据(3)中的方法遍历第二待分类地面点筛选出第三待分类地面点，第一至第三待分类地面点统称待分类地面点。
构建曲面格网金字塔R1、R2、R3的具体步骤如下：
首先，以初始的gw以及t并利用(3)中过程计算得到待分类地面点集G；对其利用局部最小值方法计算更新后的地面种子点集(该步骤目的为了防止线性方程组计算病态问题)，根据更新后的点集建立新的曲面格网，按照公式(13)的曲率约束条件得到新的待分类地面点集Gnew，该过程重复3次结束该级曲面格网R1构建；然后，令gw＝gw/2且t＝t+0.1，构建下一层级的曲面格网R2，对上一步得到的待分类地面点集再次利用局部最小值方法计算更新后的地面种子点集Gnew，根据更新后的点集建立新的曲面格网，按照公式(13)的曲率约束条件得到新的待分类地面点集，该过程重复3次结束该级曲面格网R2构建；最后，令gw＝gw/2 且t＝t+0.1，构建下一层级的曲面格网R3，对上一步得到的待分类地面点集再次利用局部最小值方法计算更新后的地面种子点集Gnew，根据更新后的点集建立新的曲面格网，按照公式 (13)的曲率约束条件得到新的待分类地面点集，该过程重复3次结束该级曲面格网R3构建。
第三步：相对高程计算与高斯混合模型分类。
(1)遍历所有的待分类点，计算每一个待分类点的相对高程，将ΔH代入二元高斯混合模型(GMM)中，计算高斯混合模型对每个待分类点的相对高程的响应度，应用期望最大化算法(EM)迭代更新计算各分模型参数值αl,μl,σl 2,l＝1,2；其中，第k个待分类点qk的相对高程ΔHk＝H–Hgrid，H为待分类点的高程值，Hgrid为待分类点所在R3的曲面格网单元的高程；
相对高程的计算是为了更显著地区分地面点与非地面点。因为，一般的算法都假设如果一个点高差大于其邻居点的高差，则将其分为非地面点。但是，这一假设当遇到复杂的地形会失效如陡崖等。通过计算每个待处理点qk的相对高程ΔHk，可以有效解决这一问题，其公式如下：
ΔHk＝H–Hgrid (14)
式中，H是待处理点qk在z方向的坐标值，Hgrid是是待处理点qk所对应的R3(m,n)的曲面高程值。
将ΔHk代入二元高斯混合模型(GMM)中，可以表示每个待处理点qk的条件概率：
式中，αl为模型的混合系数且大于等于0；ΔHk是待处理点qk的相对高程；φ()是高斯密度函数，θl＝(μl,σl)是其参数。μl是待处理点qk相对高程的平均值；σl是待处理点qk相对高程的标准差；其中高斯密度函数为：
应用式(15)计算地面点的条件概率P(qk|G)，仅用到一组参数(α1,μ1,σ1)：
P(qk|G)＝α1φ(ΔHk|μ1,σ1) (17)
应用式(15)计算非地面点的条件概率P(qk|NG)时用另一组参数(α2,μ2,σ2)：
P(qk|NG)＝α2φ(ΔHk|μ2,σ2) (18)
由于高斯混合模型中既存在观测变量ΔHk，也存在隐变量G和NG，所以常规的最大似然方法无法被直接用来估计参数。此时，期望最大化算法(EM)算法提供了一个求解隐变量近似解的方案。通常，期望最大化算法包含四个主要步骤：
i.初始化输入参数，包括αl,μl,σl,l＝1,2.
首先，我们根据多级曲率约束得到的地面种子点数量以及总待处理点数量Nc，初始化α1和α2。需要注意的是，α1+α2＝1。而μ1,μ2,σ1,σ2则分别由地面种子点与非地面种子点的相对高程计算。
式中，ΔHk的含义与公式(15)一致；Nc是待处理点云总数量。
ii.E步：计算高斯混合模型对每个点的隶属度
可以分别表示为P(G|qk)和P(NG|qk)。
iii.M步：计算新一轮迭代中高斯混合模型中的参数(αl,μl,σl 2,l＝1,2)。
迭代中的M步是为了寻找使得E步中隶属度期望最大化的参数值，αl,μl,σl 2的计算公式如下：
式中，ΔHk的含义与公式(15)一致。
iv.判断期望最大化算法的收敛条件。
期望最大化算法(EM)重复(i)-(iv)步骤，直到模型参数(αl,μl,σl 2,l＝1,2)不再变化为止。
(2)由于假设地面点、非地面点的相对高程均符合高斯分布，所以根据计算出的各分模型高斯参数值，得到包含隐变量的条件概率估计值，进而使用贝叶斯公式及全概率公式求解用于区分地面点的后验概率P(G|qk)，若后验概率P(G|qk)大于等于0.5，则认为第k个待分类点qk属于地面点，并以之参与数字高程模型DEM构建；若后验概率P(G|qk)小于0.5，则认为qk属于非地面点，并舍弃该点；
由期望最大化算法(EM)计算得到模型参数后，首先，直接利用高斯混合模型公式(15-18) 计算地面点的条件概率的估计P(qk|G)，及非地面点的条件概率估计P(qk|NG)；然后，根据贝叶斯公式，待处理点qk的属于地面点的后验概率P(G|qk)及属于非地面点的后验概率P(NG |qk)的计算公式如下：
式中，P(G)代表处理点qk属于地面点的先验概率；P(NG)代表处理点qk属于非地面点的先验概率；P(qk)由全概率公式计算，P(qk)＝P(qk|G)P(G)+P(qk|NG)P(NG)。
根据条件P(G|qk)≥0.5，若满足则将待分类点qk为标记为地面点；否则，则将其分为非地面点，并舍弃之。
因此，本发明的基于多级曲率约束的高斯混合模型点云滤波方法是一种基于概率论与数理统计的分类方法，在后验概率的指导下，完成机载LiDAR点云的滤波。附图4中的(a) 至(c)展示了一个ISPRS数据集实测点云样本的过滤结果，其中，(a)为是使用原始LiDAR 点云构建的数字表面模型(DSM)，(b)是使用本发明的基于多级曲率约束的高斯混合模型点云滤波方法得到的数字高程模型(DEM)；(c)是使用经典数学形态学滤波方法得到的数字高程模型(DEM)。
从结果可以看出，经典数学形态学滤波方法得到的数字高程模型(DEM)在地形坡度变化剧烈且混合有低矮植被的区域滤波效果不佳，导致许多遗漏误差，一些地面点测量被误分类为非地面点导致错误的DEM结果。而本发明的基于多级曲率约束的高斯混合模型点云滤波方法不仅给出了一组较通用的算法参数值，且滤波精度有着明显的改善。
以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。
Claims (8)
1.基于多级曲率约束的高斯混合模型点云滤波方法，其特征在于，包括以下步骤：
步骤一，输入机载LiDAR点云，根据大数定律，首先计算每个点云到其K近邻点的平均距离，然后基于大于M倍标准差为粗差的原则，剔除不满足条件的点云，得到待处理点云；
步骤二，对待处理点云所在区域建立网格索引，并由局部最小值方法获得的地面种子点，利用薄板样条函数TPS进行插值，剔除与插值曲面高差大于曲率阈值的点，得到待分类地面点；
步骤三，计算每一个待分类地面点的相对高差，应用二元高斯混合模型GMM进行地物与地面的分类，使用期望最大化算法EM得到模型参数的最优解，并根据计算得到后验概率判断该点是否属于地物完成点云滤波。
2.根据权利要求1所述的基于多级曲率约束的高斯混合模型点云滤波方法，其特征在于，所述步骤一的具体过程为：
(1)输入机载LiDAR点云，并根据研究区域点云分布，由用户给定离群点去噪算法SOR参数，SOR参数包括近邻点个数K以及标准差乘系数M；
(2)根据大数定律，遍历每个点云计算到其K近邻点的平均距离，再基于大于M倍标准差为粗差的原则，剔除不满足
的点云，得到待处理点云；其中，
为第i个点云到其K近邻点的平均距离，KNN平均距离
标准差
n为输入的机载LiDAR点云数。
3.根据权利要求1所述的基于多级曲率约束的高斯混合模型点云滤波方法，其特征在于，所述步骤二的具体过程为：
(1)先对步骤一得到的待处理点云，由局部最小值法获取初始地面种子点集Gseed；给定格网大小gw，构建网格索引，根据网格行列号(m,n)计算出网格单元中心的真实平面坐标值(x,y)，再由最小化薄板样条能量函数Ef的优化原则，选择待插值格网点的N个近邻地面种子点作为控制点；
(2)令控制点权重满足插值函数f(x,y)存在二阶偏导数约束条件，根据拉格朗日乘数法，解算出插值函数f(x,y)中的未知系数，进而插值出网格单元中心高程值；
(3)给定曲率阈值参数t，由gw和t构建曲面格网金字塔R1，并使用一个3×3的均值移动窗口经过每个待处理点云所在格网，得到一个平滑后的曲面格网Rnew，根据Rnew高程值znew(m,n)进而计算平均曲率c(m,n)＝znew(m,n)+t；遍历待处理点云，若第j个待处理点云Pj(xj,yj,zj)的高程值zj大于等于c(m,n)，则将其加入非地面点，否则将其作为第一待分类地面点进入(4)；
(4)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔R2，并根据(3)中的方法遍历第一待分类地面点筛选出第二待分类地面点，进入(5)；
(5)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔R3，并根据(3)中的方法遍历第二待分类地面点筛选出第三待分类地面点，第一至第三待分类地面点统称待分类地面点。
4.根据权利要求3所述的基于多级曲率约束的高斯混合模型点云滤波方法，其特征在于，(1)中根据下式换算真实平面坐标值(x，y)：
式中，(xmin，ymin)为输入的机载LiDAR点云的最小边界盒的x、y方向上最小坐标值。
5.根据权利要求3所述的基于多级曲率约束的高斯混合模型点云滤波方法，其特征在于，(2)中插值函数：
式中，a1、a2、a3为反映插值曲面总体趋势的系数；ws表示与第s个控制点cs相关的权重系数；
为与薄板样条函数TPS相关的径向基核函数，(xs，ys)为cs的x、y方向坐标值，N为控制点个数；
f(x,y)的二阶偏导数约束为：
式中，zs为cs在z方向的坐标值；
根据拉格朗日乘数法，利用以下线性方程组求解f(x,y)中的未知系数：
式中，矩阵U的第I行第J列元素
(xI，yI)和(xJ，yJ)分别为第I控制点cI和第J个控制点cJ的x、y方向坐标值；
最后，待插值单元高程值z的计算公式为：z＝f(x,y)。
6.根据权利要求1所述的基于多级曲率约束的高斯混合模型点云滤波方法，其特征在于，所述步骤三的具体过程为：
(1)遍历所有的待分类点，计算每一个待分类点的相对高程并代入二元高斯混合模型GMM中，计算高斯混合模型对每个待分类点的相对高程的响应度，应用期望最大化算法EM迭代更新计算各分模型参数值αl,μl,σl 2,l＝1,2；其中，第k个待分类点qk的相对高程ΔHk＝H–Hgrid，H为待分类点的高程值，Hgrid为待分类点所在R3的曲面格网单元的高程，αl为模型的混合系数且大于等于0，μl是qk的相对高程的平均值；σl是qk的相对高程的标准差；
(2)根据计算出的各分模型高斯参数值，得到包含隐变量的条件概率估计值，进而使用贝叶斯公式及全概率公式求解qk属于地面点的后验概率P(G|qk)，若P(G|qk)大于等于0.5，则认为qk属于地面点，并以之参与数字高程模型DEM构建；若P(G|qk)小于0.5，则认为qk属于非地面点，并舍弃该点；
其中，
P(qk)＝P(qk|G)P(G)+P(qk|NG)P(NG)，P(G)是qk属于地面点的先验概率，P(NG)是qk属于非地面点的先验概率，P(qk|G)＝α1φ(ΔHk|μ1,σ1)，P(qk|NG)＝α2φ(ΔHk|μ2,σ2)，
7.根据权利要求6所述的基于多级曲率约束的高斯混合模型点云滤波方法，其特征在于，(1)中根据下式更新计算各分模型参数值αl,μl,σl 2,l＝1,2：
式中，P(NG|qk)为qk属于非地面点的后验概率；Nc是待处理点云总数量。
8.根据权利要求7所述的基于多级曲率约束的高斯混合模型点云滤波方法，其特征在于，
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011253108.3A CN112308966A (zh) | 2020-11-11 | 2020-11-11 | 基于多级曲率约束的高斯混合模型点云滤波方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011253108.3A CN112308966A (zh) | 2020-11-11 | 2020-11-11 | 基于多级曲率约束的高斯混合模型点云滤波方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112308966A true CN112308966A (zh) | 2021-02-02 |
Family
ID=74325950
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011253108.3A Pending CN112308966A (zh) | 2020-11-11 | 2020-11-11 | 基于多级曲率约束的高斯混合模型点云滤波方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112308966A (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112862844A (zh) * | 2021-02-20 | 2021-05-28 | 苏州工业园区测绘地理信息有限公司 | 基于车载点云数据的道路边界交互式提取方法 |
| CN113109793A (zh) * | 2021-03-22 | 2021-07-13 | 中交广州航道局有限公司 | 自适应分辨率水深曲面滤波方法及装置 |
| CN116197913A (zh) * | 2023-03-23 | 2023-06-02 | 广东技术师范大学 | 一种基于点云处理的机器人加工路径规划方法和存储介质 |
-
2020
- 2020-11-11 CN CN202011253108.3A patent/CN112308966A/zh active Pending
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112862844A (zh) * | 2021-02-20 | 2021-05-28 | 苏州工业园区测绘地理信息有限公司 | 基于车载点云数据的道路边界交互式提取方法 |
| CN112862844B (zh) * | 2021-02-20 | 2024-01-05 | 园测信息科技股份有限公司 | 基于车载点云数据的道路边界交互式提取方法 |
| CN113109793A (zh) * | 2021-03-22 | 2021-07-13 | 中交广州航道局有限公司 | 自适应分辨率水深曲面滤波方法及装置 |
| CN116197913A (zh) * | 2023-03-23 | 2023-06-02 | 广东技术师范大学 | 一种基于点云处理的机器人加工路径规划方法和存储介质 |
| CN116197913B (zh) * | 2023-03-23 | 2023-12-05 | 广东技术师范大学 | 一种基于点云处理的机器人加工路径规划方法和存储介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112308966A (zh) | 基于多级曲率约束的高斯混合模型点云滤波方法 | |
| CN106529469B (zh) | 基于自适应坡度的无人机载LiDAR点云滤波方法 | |
| CN106199557B (zh) | 一种机载激光雷达数据植被提取方法 | |
| CN110992341A (zh) | 一种基于分割的机载LiDAR点云建筑物提取方法 | |
| Lee et al. | Fusion of lidar and imagery for reliable building extraction | |
| CN112101430B (zh) | 用于图像目标检测处理的锚框生成方法及轻量级目标检测方法 | |
| CN111080684A (zh) | 一种点邻域尺度差异描述的点云配准方法 | |
| CN110458174B (zh) | 一种无序点云关键特征点精确提取方法 | |
| CN112101278A (zh) | 基于k近邻特征提取和深度学习的宅基地点云分类方法 | |
| CN111598780B (zh) | 一种适用于机载LiDAR点云的地形自适应插值滤波方法 | |
| CN111340723B (zh) | 一种地形自适应的机载LiDAR点云正则化薄板样条插值滤波方法 | |
| CN110335352B (zh) | 一种机载激光雷达点云的双基元多分辨率层次滤波方法 | |
| CN105118090A (zh) | 一种自适应复杂地形结构的点云滤波方法 | |
| CN115222625A (zh) | 一种基于多尺度噪声的激光雷达点云去噪方法 | |
| CN112200083B (zh) | 一种基于多元高斯混合模型的机载多光谱LiDAR数据分割方法 | |
| CN111754447A (zh) | 基于多状态上下文隐马尔科夫模型的红外和可见光图像融合方法 | |
| CN114386466B (zh) | 一种用于脉冲星搜寻中候选体信号挖掘的并行的混合聚类方法 | |
| CN110111300B (zh) | 一种图像变化检测方法 | |
| CN114898118A (zh) | 基于多源点云的输电线路房屋拆迁量自动统计方法及系统 | |
| CN116186864B (zh) | 一种基于bim技术的深基坑模型快速建模方法及系统 | |
| CN111860359B (zh) | 一种基于改进随机森林算法的点云分类方法 | |
| CN109145738A (zh) | 基于加权非凸正则化和迭代重约束低秩表示的动态视频分割方法 | |
| CN108595373B (zh) | 一种无控制dem配准方法 | |
| CN115393631A (zh) | 基于贝叶斯层图卷积神经网络的高光谱图像分类方法 | |
| CN113449920A (zh) | 一种风电功率预测方法、系统及计算机可读介质 |
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 |