一种基于Snake模型多基元融合的点云滤波方法与流程

本发明涉及地理空间信息系统技术领域，特别是涉及一种基于snake模型多基元融合的点云滤波方法。

背景技术：

机载lidar(lightdetectionandranging)技术提供了一种全新的获取高时空分辨率地球空间信息的观测手段，它被誉为继全球定位系统以来在测绘遥感领域的又一场技术革命。机载lidar技术已广泛应用于道路提取、电力巡线、植被参数估测、变化检测、城市三维模型建立等生产、生活的众多领域。而实现以上点云后处理应用，一个非常关键的步骤便是点云滤波，即从lidar点云中去除地物点而保留地形点。

点云滤波是点云数据处理的基础性工作，近年来大量的国内外研究人员对点云滤波进行了针对性的研究。虽然已有大量的点云滤波算法，并且部分算法已进行商业化应用(如terrasolid软件应用不规则三角网渐进加密滤波算法)，但点云滤波算法研究仍处于发展阶段，还有许多问题没有得到解决，尤其是现有的滤波算法并未充分挖掘利用点云自身的特征信息，致使滤波算法自动化程度较低且难以适应复杂的地形环境。因此，迫切需要提出新的滤波理论与方法。

近年来，关于点云滤波算法的研究有很多，根据滤波基本处理单元(基元)类型的不同，可以将这些算法分为以下两类：点基元滤波算法和对象基元滤波算法。

其中，点基元滤波算法主要通过点与邻近点间的几何关系来判定一点是否属于地面点。根据临近点间几何关系定义的不同，点基元滤波算法具体可分为基于坡度的滤波算法、基于形态学的滤波算法和基于曲面的滤波算法。目前基于点基元的滤波算法较多，但此类型算法滤波理论相对单一，基本上都是在原有方法(基于坡度、基于形态学、基于曲面)的基础上进行改进的，缺乏新的滤波原理与方法。此外，点基元滤波算法往往可利用特征不足，并未充分挖掘点云数据本身的特征信息，使得该类型算法在复杂地形区域往往滤波效果较差。

而对象基元滤波算法通常包含两步，首先采取某种分割方法对点云进行分割，然后再对分割的结果按照某种设定的规则进行点云滤波。点云分割的方法有很多，如扫描线的分割法、区域生长法、随机抽样一致法(ransac)等。在对分割结果进行滤波判断时，大多数算法通常都基于地面点聚类区域要低于地物点聚类区域这一假设。在滤波原理相似的前提下，相较于点基元滤波算法，对象基元滤波算法往往能够获得更好的滤波效果。但此类型算法的滤波结果严重依赖点云分割的质量，这也体现了采用单一基元进行点云滤波的局限性。

虽然滤波方法有很多，但依然存在以下三个方面的问题亟待解决：

(1)大多数滤波方法无法兼顾点云的局部尺度信息与全局尺度信息，使得滤波方法在部分复杂地形区域(如陡坡、不连续地形、斜坡上建有房屋等)难以取得良好的滤波效果。

(2)采用单一基元实现点云滤波无法充分挖掘利用点云数据的自身特征，致使滤波方法鲁棒性较差，难以适应复杂地形环境。

(3)为能适应各种地形环境，大多数滤波算法往往需要进行复杂的参数设置与阈值调节，这不仅降低了算法的自动化程度，而且也不利于经验缺乏的工作人员进行滤波实现。

技术实现要素：

为此，本发的实施例提出一种基于snake模型多基元融合的点云滤波方法，解决现有技术中的上述问题。

一种基于snake模型多基元融合的点云滤波方法，包括以下步骤：

s1，加载机载雷达las格式的数据文件，获取点云xyz坐标信息；

s2，对点云数据进行格网化组织，采用形态学开运算探测点云数据中的噪声点，并进行剔除；

s3，获取初始snake模型，同时计算此时的初始能量；

s4，采用meanshift自动聚类方法，将每个点基元移动到密度函数的局部极大值点处，获取多个对象基元；

s5，利用对象基元不断更新snake模型，直至模型能量达到全局能量最小化；

s6，利用步骤s5获取的全局能量最小化的snake模型计算局部滤波阈值，并按点基元实现点云自动滤波。

根据本发明提供的基于snake模型多基元融合的点云滤波方法，采用多基元融合(点基元和对象基元)的滤波策略，可以有效避免单一基元滤波法鲁棒性差的缺陷，充分利用不同基元提供的特征信息，提高滤波方法对复杂地形区域的适应性。主要构建snake模型，通过设置具有地面点吸附力的目标能量函数，使得模型能够兼顾点云的局部尺度信息与全局尺度信息，充分挖掘点云数据本身的自动化潜力。依据该模型将滤波阈值表示为地形坡度变化的线性函数，实现滤波阈值的自适应更新，避免单一阈值滤波精度差而多阈值需要过多参数设置的缺陷，因此与现有技术相比，本发明提供的方法具有以下优点：(1)构建合理的snake模型，兼顾点云的局部尺度信息与全局尺度信息，提高滤波方法在复杂地形区域的滤波精度；(2)采用多基元融合的滤波策略，充分利用不同基元空间特性，提高滤波方法的鲁棒性；(3)基于全局能量最小化的snake模型，实现滤波阈值随地形变化的自适应更新，提升算法自动化程度。

另外，根据本发明上述的基于snake模型多基元融合的点云滤波方法，还可以具有如下附加的技术特征：

进一步地，所述步骤s2具体包括以下步骤：

s21，将点云数据进行格网组织，生成深度图像i(x,y)；

s22，对格网数据进行形态学膨胀运算，结构元素尺寸为3×3；

s23，将形态学膨胀运算前后高差变化大于阈值的孤立点判定为噪声点，并进行剔除。

进一步地，所述步骤s21具体包括以下步骤：

s211，获取点云数据的二维平面包围盒；

s212，计算平均点间距，并以平均点间距作为格网尺寸对二维平面包围盒进行格网剖分；

s213，如果格网内有多个点存在，则取格网内所有点高程值的最小值作为该格网的特征值；

s214，如果格网内没有数据，则取与该格网中心最邻近点的高程值作为该格网的特征值。

进一步地，所述步骤s3具体包括以下步骤：

s31，对去噪后的深度图像进行大尺度的形态学腐蚀运算，获取初始snake模型；

s32，计算模型的初始能量，即计算模型的内部能量与外部能量之和。

进一步地，所述步骤s32具体包括以下步骤：

s321，利用初始snake模型计算地形坡度；

s322，将模型内部能量表示成地形坡度的log函数并进行计算；

s323，计算初始模型与实际地形之间的高差；

s324，将模型外部能量表示成高差指数函数并进行计算；

s325，计算内部能量与外部能量之和获取模型的初始能量。

进一步地，所述步骤s4具体包括以下步骤：

s41，任选一个点基元，如果该点基元未被访问过，则统计在以该点基元为圆心带宽h为半径的圆形范围内的所有点基元；

s42，计算该圆形范围内所有点基元的meanshift向量mh,g(p)，如果mh,g(p)大于阈值则进行步骤s43，如果小于等于阈值则进行步骤s44；

s43，根据meanshift向量更新计算圆心坐标位置，并统计以该圆心为中心半径为h范围内的所有点基元，继续执行步骤s42；

s44，结束迭代循环，记录该点基元对应的模态点位置，并执行步骤s41，直到所有的点基元都被访问过位置；

s45，将具有相同或相近模态点的点基元标记为同一类，获取多个对象基元。

进一步地，所述步骤s5具体包括以下步骤：

s51，获取步骤s4中的对象基元，并依次更新snake模型，计算更新后的模型能量e′(m)；

s52，计算snake模型更新前后的能量差，并按metropolis准则接受模型更新；

s53，重复步骤s51-s52，直到不再有新的对象基元可以使模型能量降低为止。

其中，所述步骤s51的具体方式如下：获取任一步骤s4中的对象基元，如果该对象基元没有被访问过，则利用该对象基元对snake模型相应位置进行替换更新，并将该对象基元标记为已访问。

进一步地，所述步骤s52具体包括以下步骤：

s521，初始化温度常量t以及常数k；

s522，计算snake模型更新前后的能量差δe；

s523，如果δe小于0，则接受更新后的模型作为当前的snake模型。如果δe大于等于0，则按照metropolis准则以概率p(δe)＝exp(δe/(kt))接受该模型；

s524，常量温度减少t，且t→0。

进一步地，所述步骤s6具体包括以下步骤：

s61，依次遍历各个点基元，获取各点基元的二维格网坐标(i,j)；

s62，依据步骤s5所获取的snake模型以及各点基元所对应的二维格网坐标获取各点基元所对应的滤波阈值δh；

s63，计算各点基元到snake模型的距离ndsm；

s64，如果ndsm小于等于δh，则将该点基元标记为0，即地面点，如果ndsm大于δh，则将该点基元标记为1，即地物点。

进一步地，所述步骤s62具体包括以下步骤：

s621，利用snake模型计算地形的局部梯度；

s622，将滤波阈值表示为梯度变化的线性函数；

s623，依据各基元所对应的二维格网坐标(i,j)获取各点基元所对应的滤波阈值。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实施例了解到。

附图说明

本发明实施例的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1是根据本发明一实施例的基于snake模型多基元融合的点云滤波方法的流程图；

图2是ⅰ预处理阶段的流程图；

图3是ⅱ对象基元更新模型阶段的流程图；

图4是ⅲ全局能量最小化处理阶段的流程图；

图5是ⅳ点云滤波阶段的流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1，本发明一实施例提出的基于snake模型多基元融合的点云滤波方法，主要深入挖掘利用不同基元提供的特征信息，提出多基元融合的滤波策略，提高滤波方法对复杂地形区域的适应性以及滤波精度。所构建的snake模型能够兼顾点云的局部尺度信息和全局尺度信息，能够充分发掘点云数据本身的自动化潜力。在此基础上，利用snake模型实现滤波阈值的自适应更新，避免单一阈值滤波精度差而多阈值需要过多参数设置的缺陷。

本实施例中，涉及的方法实施步骤流程图如图1所示，具体包括以下步骤：

s1，加载机载雷达las格式的数据文件，获取点云xyz坐标信息；

s2，对点云数据进行格网化组织，采用形态学开运算探测点云数据中的噪声点，并进行剔除；

s3，获取初始snake模型，同时计算此时的初始能量；

s4，采用meanshift自动聚类方法，将每个点基元移动到密度函数的局部极大值点处，获取多个对象基元；

s5，利用对象基元不断更新snake模型，直至模型能量达到全局能量最小化；

s6，利用步骤s5获取的全局能量最小化的snake模型计算局部滤波阈值，并按点基元实现点云自动滤波。

其中，上述六个步骤构成四个阶段：ⅰ预处理阶段；ⅱ对象基元更新模型阶段；ⅲ全局能量最小化处理阶段；ⅳ点云滤波阶段。下面就这四个阶段进行具体说明。

ⅰ预处理阶段包括图1中的步骤s1～s3，具体的，预处理阶段算法流程图如图2所示，具体包括以下步骤：

加载机载雷达las格式的数据文件，获取点云xyz坐标信息；

对点云数据进行格网化组织，采用形态学开运算探测点云数据中的噪声点，并进行剔除；

获取初始snake模型m(x,y)，同时计算此时的初始能量e(m)。

其中，步骤s2采取以下方法进行具体实现：

首先，将格网边长设定为点云平均点间距。将格网边长设定为平均点间距的目的在于确保大部分格网内都有点存在，减少空白格网的数目，尽可能小地减小因对空白格网内插而引入的误差。

继而采用下式计算各个点的格网坐标：

式中，(xk,yk)为第k个点的xy坐标，(xmin,ymin)为点云xy方向的最小值。如果有多于一个点落入同一个格网，那么该格网的特征值取这些点高程值的最小值。点云格网化组织完成后，采用形态学开运算来探测点云噪声点并进行剔除，具体实现方式如下：

首先，对格网化的二维点云数据进行结构元素为3×3形态学膨胀运算。形态学膨胀运算即为取结构元素内高程的最大值最为该格网新的特征值，公式表示如下：

zd(i,j)＝max{zori(i+u,j+v)|(u,v)∈3×3}

接下来，依次遍历点云中各个点，找到各个点所对应的格网，计算各个点与该点所对应格网膨胀后的高程值之间的差值，将差值大于阈值的孤立点判定为噪声点，并进行剔除。

式中，t1为形态学膨胀运算前后高差变化阈值，可设为5m；t2为与该点相邻点高程阈值，可设为1m；t3为与该点向临近点个数阈值，可设为3个。上述公式表明，将形态学膨胀运算前后高程大于5m的孤立点(该点的临近点个数不多于3个)判定为噪声点。

噪声剔除完成后，便可获取初始snake模型，并计算模型的初始能量。

对于三维地形环境，snake模型可以视为一个离散的二维曲面。为了便于模型能量计算，本发明将该模型以2.5维的形式进行表示，定义为m(i,j)。其中，(i,j)表示模型对应格网化点云的横纵坐标值，m(i,j)为格网对应的特征值。初始采用大尺度形态学腐蚀运算获得，即m(i,j)＝min{zori(i+u,j+v)|(u,v)∈30×30}。形态学腐蚀运算即使取结构元素内高程值的最小值作为该格网新的特征值：

snake模型主要由目标能量函数控制，目标能量函数又由内部能量函数(内部力)和外部能量函数(外部力)组成。为充分挖掘利用点云局部特征和全局特征，需设置科学合理的目标能量函数。本发明设置的目标能量函数表示如下：

e(m)＝∑(eint(m)+eext(m))

式中，eint(m)表示内部能量函数，eext(m)表示外部能量函数。内部能量函数控制模型的平滑性和连续性；外部能量函数则约束模型与全局地形最佳贴合。本发明将内部能量函数表示为地形坡度的函数，公式表示如下：

式中，g表示坡度，gsd为格网尺寸。由上式可以看出，内部能量函数与模型坡度变化成正相关，当snake模型m(i,j)较为平滑、起伏不大时，eint(m)取值较小。当snake模型为一平面时，eint(m)取到极小值。

外部能量函数产生外部力，旨在使snake模型m(i,j)能够吸附于实际地形i(i,j)。本发明将外部能量函数表示为模型与实际地形间高差的函数，公式表示如下：

式中，i(i,j)表示实际地形，δ表示高差。由公式(3)可以看出，当snake模型m(i,j)越靠近实际地形i(i,j)时，外部能量eext(m)越小。当m(i,j)≡i(i,j)时，eext(m)取到极小值。

综合内部能量公式和外部能量公式可以发现，当snake模型靠近实际地形时，虽然外部能量eext(m)会变小，但与此同时，snake模型本身会变得不那么平滑，模型坡度会变化较大，相应的内部能量eint(m)会变大。因此，目标能量函数e(m)最小化的过程即是不断调和平衡内部能量和外部能量的过程。当e(m)达到最小值时，此时对应的snake模型既能保证模型本身相对较为平滑，又能保证模型能够吸附于全局地面点，公式表示如下：

m(i,j)＝argmin(e(m))

ⅱ对象基元更新模型阶段，包括图1中的步骤s4，具体的，算法流程如图3所示，具体步骤如下：

利用meanshift方法能够将具有相近属性的点云聚集在一起，整体表现为多个对象基元。继而将这些对象基元均标记为未访问，公式表示如下：

o^k(i,j)＝0,k＝1,2,…,n

式中，o^k(i,j)表示第k个对象基元，n为meanshift分割后对象基元的总数。

选取一个对象基元，如果该对象基元未被访问过，则用该对象基元与snake模型相应位置进行更新，生成新的snake模型m′(i,j)，公式表示如下：

接着，根据新模型m′(i,j)计算模型能量e′(m)。

ⅲ全局能量最小化处理阶段包括图1中的步骤5，具体的，算法流程如图4所示，具体步骤如下：

本发明赋予对象基元一种时变且最终趋于零的概率突跳性，来获取全局能量最小值。它能够以一定的概率来接受一个比当前解要差的解，因此有可能会跳出局部的最优解，达到全局的最优解。接受概率主要基于metropolis准则计算得到，该准则通常表示如下：

metropolis准则表明，当更新后的模型能量变大时，出现能量差为δe的对象基元接受概率为p(δe)＝exp(δe/(kt))，其中为常数k常数，exp表示自然指数。具体能量最小化过程表述如下：

①任选一个未被访问过的对象基元o^k(i,j)，获取新模型m′(x,y)，并计算新模型对应能量e′(m)。

②计算模型m′(i,j)和m(i,j)间的能量差δe，如果δe小于0，则接受该对象基元；如果δe不小于0，则按metropolis准则接受该对象基元。

③如果接受对象基元，则用该对象基元更新snake模型m(i,j)，并将该对象基元标记为已访问，即o^k(i,j)＝1。

④重复执行步骤①-③，直到不再有新的对象基元可以使模型能量降低为止，即达到了模型全局能量最小化的状态，获取最终的snake模型m(i,j)。

ⅳ点云滤波阶段包括图1中的步骤s6，具体的，算法流程如图5所示，具体步骤如下

采用点基元进行滤波判断可以有效地避免滤波精度对点云分割结果的过度依赖最终获取的snake模型能够反映全局地形的变化趋势，也能够体现地形局部起伏特征。为实现算法自动化滤波，需借助该模型计算出局部梯度反应地形起伏变化。本发明采用下公式计算局部梯度：

局部梯度计算出来后，可将滤波阈值表示为局部梯度的线性函数，即δh＝f(gradient)，从而实现滤波阈值随地形变化而自适应调节。继而，计算各点基元到snake模型的距离，本发明将其定义为正则化高程，用ndsm来表示。如果ndsm小于滤波阈值δh，则该点基元为地面点；否则，该点基元判定为非地面点并进行剔除。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。