CN115294302A - 一种基于断裂线约束的机载点云快速滤波方法 - Google Patents

Abstract

本发明公开了一种基于断裂线约束的机载点云快速滤波方法，包括以下步骤：(1)采集测量区域内机载激光点云数据，根据各点云及其邻域点平均曲率值、高程值提取断裂点；(2)获取包含点云数据的扫描线；(3)针对每条扫描线进行粗过滤以滤除扫描线上属于非地面点与地物点的点云并提取扫描中心和基准点；所述扫描中心点为粗过滤后扫描线最左端的点云，基准点为距离扫描中心点最近的点云；(4)针对每条扫描线进行精细过滤；(5)采用二次曲面拟合法滤除测量区域内属于非地面点的点云，获取最终的滤波后的点云。本发明构建了考虑断裂线约束的扫描线滤波模型，适用于多种复杂场景并考虑高程突变、噪声点的影响，提高地面点提取的完整度和准确度。

Description

一种基于断裂线约束的机载点云快速滤波方法

技术领域

本发明涉及激光点云测量数据智能处理与地理信息数据应用研究的技术领域，具体是涉及一种基于断裂线约束的机载点云快速滤波方法。

背景技术

机载激光点云滤波作为准确分离和识别地面点的手段，是建立高精度数字地面模型的关键环节。准确、高精度的点云滤波对地面信息提取、道路规划以及目标识别精度的提高等应用具有重要作用。激光扫描系统能够通过扫描测得方位信息，为大面积测绘提供了一种快速准确的方案。在点云数据处理过程中，地面信息十分重要，但地面点数量巨大，其完整性易受噪声、遮挡等因素影响，这些因素造成数据处理速度低，目标分类和识别困难。

目前，基于机载激光点云滤波方法主要有三类：第一类是基于坡度变化，逐点计算点与相邻点之间的坡度，根据两点间实际坡度与坡度阈值比较滤除地物点。这类方法在处理陡坡、高程突变区域时，会错误滤除大量地面点，且坡度越陡，过滤精度越低。第二类是数学形态学方法，对地形点云进行腐蚀膨胀，构造开运算，利用开运算后非地面点高程变化大这一特点，将高差大于阈值的点分类为非地面点并进行滤除。这类方法得到的滤波结果依赖于最佳窗口尺寸的选择，较小的窗户尺寸会保留较大的建筑物，较大的窗口会平滑地形细节。第三类基于曲面的方法，通过迭代地从原始数据集中选择地面测量值，逐渐逼近地面，创建一个近似裸露地球的表面。此类方法在地形变化平缓的地区效果比较好，但对于大型建筑物往往会误把其顶面错判为地形表面。

近年来，基于激光雷达扫描线理论广泛应用于机载激光点云滤波，基于扫描线的滤波模型能顾及整个扫描区域地形的特征，能够快速高效处理大面积、大数据量区域，为机载激光点云中地物点的滤除提供了一种新思路。Yin等用球面坐标代替笛卡尔坐标，算法不假设道路位于最低位置，利用方位角和径向距离球面坐标作为判断标准，通过球坐标中径向距离曲线的断点和转折点分割地面点和非地面点，该算法可以避免过度分割和欠分割。Zheng等基于可变半径圆环和B样条拟合处理点云，得到地面点云。针对每条扫描线，利用半径可变的圆环滚过，保留扫描线上被圆环滚过的点，采用B样条拟合曲面，遍历区域中每个点，根据拟合高程与实际高程的差值实现保留地面点并滤除地物点的目的。JM Sánchez等提出一种基于扫描线处理的地面滤波算法，处理核心是利用一维迭代样条插值算法来改进地面点的提取。

但是上述提及的基于一维地形特征的滤波模型多适用于地形平坦时地面点的提取，难以适用于地面崎岖、场景复杂的区域，不能避免由断裂线等引起高程突变从而导致大量地面点被错误滤除的情形。

发明内容

发明目的：针对以上缺点，本发明提供一种基于断裂线约束的机载点云快速滤波方法，能够基于扫描线滤波的优势，构建了考虑断裂线约束的扫描线滤波模型，适用于多种复杂场景并考虑高程突变、噪声点的影响，提高地面点提取的完整度和准确度。

技术方案：为解决上述问题，本发明公开一种基于断裂线约束的机载点云快速滤波方法，具体包括以下步骤：

(1)采集测量区域内机载激光点云数据，根据各点云与其对应邻域点云的平均曲率值、高程值提取断裂线上的断裂点；

(2)根据测量区域的范围和预设的扫描线间距确定各条扫描线的横坐标；获取各点云的横坐标与纵坐标，计算各点云在横轴方向上到各条扫描线的距离，将位于设定距离内的点云归入对应的扫描线，再将扫描线上各点云按照纵坐标值升序排列；

(3)针对每条扫描线进行粗过滤以滤除扫描线上属于非地面点与地物点的点云并提取扫描中心和基准点；所述扫描中心点为粗过滤后扫描线最左端的点云，基准点为距离扫描中心点最近的点云；

(4)针对每条扫描线进行精细过滤，所述精细过滤具体包括：

(4.1)按照扫描线上的点云顺序获取基准点B相邻的点云为待测点，判断该待测点是否为断裂点，若为断裂点则保留该点云；若为非断裂点，则以该待测点构建有向向量并计算向量夹角ρ，若ρ不大于阈值，则保留该待测点否则滤除该待测点；具体公式为：

式中，B为基准点,P_i为待测点；O为扫描中心点；

(4.2)若上述待测点被保留，则将该保留的待测点替换为基准点，重复步骤(4.1)直至遍历扫描线上所有的点云；若上述待测点被滤除，则不替换该基准点B，重复步骤(4.1)直至遍历扫描线上所有的点云；

(5)采用二次曲面拟合法滤除测量区域内属于非地面点的点云，获取最终的滤波后的点云。

进一步的，步骤(1)中所述提取断裂点具体包括以下步骤：

(1.1)计算测量区域中每个点云与其邻域内点云的平均曲率值H；具体公式为：

式中，E＝r_x·r_x、F＝r_x·r_y、G＝r_y·r_y、

r_y，r_xx，r_xy，r_yy均为z(x,y)的偏导数；z(x,y)＝ax²+bxy+cy²为任意一点云及其邻域内的点云使用最小二乘法拟合曲面并对其进行参数化后的表达式，a、b、c均为拟合参数；

(1.2)针对测量区域中每个点云，计算高程小于该点的邻域点数与全部领域点数的比值G；

(1.3)根据各点云的H和G计算各点云的特征值D，保留特征值D大于预设值的点云作为断裂点；具体公式为：

D＝|H|q₁+Gq₂

式中，q₁为Ⅰ类断裂点权重，q₂为Ⅱ类断裂点权重。

进一步的，步骤(3)中所述的粗过滤具体包括以下步骤：

(3.1)定义每条扫描线上最左端的点云为初始点P₁，按照扫描线上的点云顺序选取与初始点相邻的点云，计算该点云与初始点的高度差，若高度差小于设定高度差阈值，则该点云为地面点记为P₂；否则滤除该点继续按照扫描线上的点云顺序选取后续新点重复上述操作，直至确定地面点P₂；

(3.2)选取与P_i+1相邻的点云记为P_i+2，以P_i P_i+1、P_i+1P_i+2长度分别构建等边三角形P_i O_i P_i+1、等边三角形P_i+1O_i+1P_i+2；若∠O_iP_i+1O_i+1小于夹角阈值，则判定P_i+2点为地面点；若∠O_iP_i+1O_i+1大于夹角阈值，则判定该P_i+2点为非地面点，滤除该P_i+2点并继续按照扫描线上的点云顺序选择后续新点作为新的P_i+2，重新构建等边三角形并获取夹角∠O_iP_i+1O_i+1，再将∠O_iP_i+1O_i+1与夹角阈值比较直至获取为地面点的P_i+2点；i取1……n；

(3.3)令i自加1，重复步骤(3.2)直至后续新点为扫描线上最后一个点云，滤除扫描线上的非地面点，确定扫描线上所有地面点；

(3.4)计算扫描线上每个点云的点云平坦度，保留点云平坦度大于平坦度阈值的点云。

进一步的，步骤(5)具体包括以下步骤：

(5.1)根据原始点云数据中每个点云的横纵坐标范围划分格网，不同网格之间形成重叠区域、非重叠区域；获取落入单元格内的滤波后扫描线上的m个地面点，采用最小二乘原理求得误差平方和σ²达到最小时的单元格的曲面参数A′、B′、C′、D′、E′、F′；具体的公式为：

式中，z(x_i,y_i)为点云的坐标；

(5.2)确定每个格网曲面的参数后，判断原始点云中每个点云所在的格网位置：若原始点云落入格网非重叠区域，则计算该点云在单个格网的拟合高程z(x,y)；若原始点云落入格网重叠区域，则计算该点云在不同格网曲面拟合高程z(x,y)的平均值；z(x,y)获取的公式为：

z(x,y)＝A'x²+B'y²+C'xy+D'x+E'y+F'

(5.3)比较各点云的拟合高程值与实际高程的差值，滤除差值大于预设阈值的点云，获得最终的滤波后的点云。

进一步的，步骤(3)中所述的粗过滤还包括：

(3.5)对于点云中的任一点云，获取高程大于该点的邻域点数与全部邻域点数的比值，保留比值小于设定阈值的点云。

进一步的，步骤(3.2)中由等边三角形几何关系计算获得O_i的坐标定义具体公式为：

式中，P_i表示为

P_i+1表示为

P_i+2表示为(x_i+2 ^p，y_i+2 ^p，z_i+2 ^p)，O₁表示为

进一步的，步骤(3.4)中点云平坦度的计算公式为：

式中，X_i表示为点云的领域，0＜i≤n；D表示为协方差矩阵；

为空间邻域的中心点；λ为点云平坦度；λ₁、λ₂、λ₃分别为协方差矩阵D在三维空间的特征值。

此外，本发明还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述任一方法的步骤。一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一方法的步骤。

有益效果：本发明所述一种基于断裂线约束的机载点云快速滤波方法相对于现有技术，其显著优点是：联合曲率值和高程差获取断裂点适用于不同类型测区滤波；针对各类场景的机载激光点云数据，一方面考虑构建等边三角形计算坡度值、点云平坦度进行粗过滤；另一方面加入断裂点约束，克服高程突变、噪声点的影响完成精细过滤；最后采用二次曲面拟合法滤除测量区域内的非地面点，提高地面点提取的完整度和准确度，形成数字高程模型。

附图说明

图1所示为本发明所述方法的流程图；

图2所示为本发明所述方法中不同断裂线的示意图，图2(a)为Ⅰ类断裂线的示意图，图2(b)为Ⅱ类断裂线的示意图；

图3所示为本发明所述方法中提取断裂线点的结果图；

图4所示为本发明所述方法扫描线粗过滤示意图；

图5所示为本发明所述方法中构建的向量以及向量夹角的示意图；

图6所示为本发明所述方法中经过细滤波前后的对比示意图，图6(a)为细滤波前的示意图，图6(b)为细滤波后的示意图；

图7所示为本发明所述方法中对测量区域的格网划分示例的示意图；

图8所示为本发明所述方法中不同种类复杂场景中地面点云滤波结果示意图；图8(a)为地形1滤波后的示意图，图8(b)为地形2滤波后的示意图，图8(c)为地形3滤波后的示意图，图8(d)为地形4滤波后的示意图。

具体实施方式

下面结合附图对本发明的技术方案进一步说明。

如图1所示，本发明提供的一种基于断裂线约束的机载点云快速滤波方法，具体包括以下步骤：

步骤一、采集测量区域内机载激光点云数据，结合高程差异和曲率值联合提取断裂线上的断裂点。

如图2(a)所示，断裂线包括在山脊线、山谷线等处曲面连续的Ⅰ类断裂线；如图2(b)所示，断裂线还包括曲面在人工建筑、陡坎等处不连续的Ⅱ类断裂线。由于Ⅰ类断裂线上的断裂点均为连续的，可以直接根据点云的平均曲率的绝对值提取Ⅰ类断裂点。平均曲率的绝对值反映了地表的起伏程度，当该值较大时，表明局部地表的起伏程度较大，反之则表明地面较平坦，能够通过提取平均曲率的绝对值较大的点作为潜在断裂点。由于在Ⅱ类断裂线处的地面点不连续且存在显著的高程差异，直接根据曲率提取该类断裂线效果欠佳，可通过统计高程小于该点的邻域点与全部领域点数的比值提取Ⅱ类断裂点，若比值大于给定阈值，则判断该点为潜在断裂点。

因为对于不同类型的断裂线上的断裂点具有不同的提取方式，为了适用于包含不同类型断裂线的测量区域的滤波，本方法联合高程差异和曲率值联合提取断裂点，具体包括以下步骤：

(1)针对任意一点云及其邻域内的点云使用最小二乘法拟合曲面并对其进行参数化获得：z(x,y)＝ax²+bxy+cy²，a、b、c均为拟合参数；z(x,y)的偏导数分别为r_x、r_y、r_xx、r_xy、r_yy；计算测量区域中每个点云与其邻域内点云的平均曲率值H；具体公式为：

式中，E＝r_x·r_x、F＝r_x·r_y、G＝r_y·r_y、

L、M、N为曲面第二基本不变量；E、F、G为曲面第二基本不变量；

为曲面主法线方向；

(2)针对测量区域中每个点云，计算高程小于该点的邻域点数与全部领域点数的比值G；

(3)根据各点云的H和G计算各点云的特征值D，公式为：

D＝|H|q₁+Gq₂

式中，q₁为Ⅰ类断裂点权重，q₂为Ⅱ类断裂点权重。

将各点云的特征值D与预设值作比较，若特征值D大于预设值，本实施例中预设值设定为0.6，则该点云作为断裂点，否者为非断裂点。根据以上步骤最终获得各断裂点形成的断裂线如图3所示，且结果显示断裂线基本符合实际情况，可满足滤波约束要求。

步骤二、对于全区域点云进行一维地形特征点提取，获取内各点云的横坐标与纵坐标。根据测量区域的范围和预设的扫描线间距确定各条扫描线的横坐标；计算各点云在横轴方向上到各条扫描线的距离，将位于设定距离内的点云归入对应的扫描线，同时将扫描线上各点云按照纵坐标值升序排列。

步骤三、对于每条扫描线进行粗过滤从而滤除扫描线上属于非地面点与地物点的点云，选取扫描中心。所述的粗过滤包括很多种技术方案，例如：通过判断相邻点的高差和坡度是否均在预设的阈值范围内进而判断相邻点是否为非地面点。本方法所述的粗过滤为先利用扫描线上的点顺序构建等边三角形计算坡度过滤多数地物点，再根据地面平坦度剔除潜在建筑物点，最后根据高程差异过滤非地面点。具体步骤包括：

(1)利用扫描线上的点顺序构建等边三角形计算坡度过滤多数地物点。

(1.1)定义每条扫描线上最左端的点云为初始点P₁，按照扫描线上的点云顺序选取与初始点相邻的点云，计算该点云与初始点的高度差，若高度差小于设定高度差阈值h，则该点云为地面点记为P₂；否则滤除该点继续按照扫描线上的点云顺序选取后续新点重复上述操作，直至确定地面点P₂；

(1.2)选取与P_i+1相邻的点云记为P_i+2，以P_i P_i+1、P_i+1 P_i+2长度分别构建等边三角形P_i O_i P_i+1、等边三角形P_i+1 O_i+1 P_i+2；若∠O_iP_i+1O_i+1小于夹角阈值，则判定P_i+2点为地面点；若∠O_iP_i+1O_i+1大于夹角阈值，则判定该P_i+2点为非地面点，滤除该P_i+2点并继续按照扫描线上的点云顺序选择后续新点作为新的P_i+2，重新构建等边三角形并获取夹角∠O_iP_i+1O_i+1，再将∠O_iP_i+1O_i+1与夹角阈值比较直至获取为地面点的P_i+2点；i取1……n；

具体的，如图4所示，选取与地面点P₂相邻的点云记为P₃，以P₁ P₂、P₂ P₃长度构建等边三角形P₁ O₁ P₂与等边三角形P₂ O₂ P₃；若∠O₁P₂O₂小于夹角阈值，则判定该P₃点为地面点；若∠O₁P₂O₂大于夹角阈值，则判定该P₃点为非地面点，滤除该P₃点并继续按照扫描线上的点云顺序选择后续新点作为新的P₃，重新构建等边三角形获取夹角∠O₁P₂O₂并与夹角阈值比较直至获取为地面点的P₃点；

其中，由等边三角形几何关系计算获得O₁的坐标，同理可计算出O₂坐标。具体公式为：

式中，当前点

P₃(x₃ ^p，y₃ ^p，z₃ ^p)，O₁(x₁ ^o，y₁ ^o，z₁ ^o)。

(1.3)令i自加1，重复步骤(3.2)直至后续新点为扫描线上最后一个点云，遍历扫描线上所有的点，以该方法滤除扫描线上的非地面点，确定扫描线上所有地面点。

(2)根据地面平坦度剔除潜在建筑物点。计算扫描线上每个点云的点云平坦度，保留点云平坦度大于平坦度阈值的点云。具体点云平坦度的计算公式为：

式中，X_i表示为点云的领域，0＜i≤n；D表示为协方差矩阵；

为空间邻域的中心点；λ为点云平坦度；λ₁、λ₂、λ₃(λ₁≤λ₂≤λ₃)分别为协方差矩阵D在三维空间的特征值。

(3)计算点云平坦度滤除扫描线上的植被、车辆、人员等地物后，需要考虑建筑物等局部邻域平坦的非地面点对扫描中心选取的影响。为剔除扫描线上的建筑物，对扫描线上点云中任一点，获取高程大于该点的邻域点数与全部邻域点数的比值，比值越小，该点为地面点的可能性越大，故保留比值小于设定阈值的点云。

针对每条扫描线进行粗过滤后，以扫描线上最左端点为扫描中心O，离扫描中心距离最近的点为基准点B。

步骤四、考虑断裂线上断裂点的影响，对于扫描线进行精细滤波。具体包括：

(1)按照扫描线上的点云顺序获取基准点B相邻的点云为待测点，判断该待测点是否为断裂点，若为断裂点则保留该点云；若为非断裂点，则以该待测点构建有向向量并计算向量夹角ρ，若ρ不大于阈值，则保留该待测点否则滤除该待测点；具体公式为：

式中，B为基准点,P_i为待测点；O为扫描中心点；

(2)若上述待测点被保留，则将该保留的待测点替换为基准点，重复步骤(1)直至遍历扫描线上所有的点云；若上述待测点被滤除，则不替换该基准点B，重复步骤(1)直至遍历扫描线上所有的点云。

如图5所示，当前基准点为P₅点，构建有向向量

获取向量夹角

该向量夹角小于阈值，认定为非地面点滤除，即P_i-1被认定为非地面点；保持该基准点不换，继续选取后续新点，重新构建有向向量，继续判定新待测点是否为地面点，直至完成扫描线上所有点云的判定。如图6所示，图6(a)为原始的点云数据图像，图6(b)为经过精细滤波后的点云图像数据，可以区分图像中的地面点和非地面点。

步骤五、采用二次曲面拟合法滤除测量区域内属于非地面点的点云，即对于扫描线之间的点云进行筛选过滤，获取最终的滤波后的点云，形成数字高程模型。

(1)根据原始点云数据中每个点云的预设的横纵坐标范围划分格网，格网单元尺寸原则上应大于最大地物的边长。不同网格之间形成重叠区域、非重叠区域；如图7(a)所示，测量区域内所有网格划分后存在非重叠区域、二格重叠区域、四格重叠区域。

获取落入单元格内的滤波后扫描线上的m个地面点，采用最小二乘原理求得误差平方和σ²达到最小时的单元格的近似曲面参数A′、B′、C′、D′、E′、F′；具体的公式为：

式中，z(x_i,y_i)为点云的坐标；

(2)确定每个格网曲面的参数后，判断原始点云中每个点云所在的格网位置：若原始点云落入格网非重叠区域，则计算该点云在单个格网的拟合高程z(x,y)；若原始点云落入格网重叠区域，则计算该点云在不同格网曲面拟合高程z(x,y)的平均值；z(x,y)获取的公式为：

z(x,y)＝A'x²+B'y²+C'xy+D'x+E'y+F'

(3)比较各点云的拟合高程值与实际高程的差值，滤除差值大于预设阈值的点云，获得最终的滤波后属于地面点的点云。

如图8所示，采用本发明所述方法针对不同的类型的测量区域进行地面点提取；如图8(a)、8(b)、图8(c)、8(d)分别为四种不同地形的滤波后的结果图，根据图中显示可以提取观测区域内所有属于地面点的点云数据。

Claims (9)

1.一种基于断裂线约束的机载点云快速滤波方法，其特征在于，包括以下步骤：

(1)采集测量区域内机载激光点云数据，根据各点云与其对应邻域点云的平均曲率值、高程值提取断裂线上的断裂点；

(2)根据测量区域的范围和预设的扫描线间距确定各条扫描线的横坐标；获取各点云的横坐标与纵坐标，计算各点云在横轴方向上到各条扫描线的距离，将位于设定距离内的点云归入对应的扫描线，再将扫描线上各点云按照纵坐标值升序排列；

(3)针对每条扫描线进行粗过滤以滤除扫描线上属于非地面点与地物点的点云并确定扫描中心和基准点；所述扫描中心点为粗过滤后扫描线最左端的点云，基准点为距离扫描中心点最近的点云；

(4)针对每条扫描线进行精细过滤，所述精细过滤具体包括：

(4.1)按照扫描线上的点云顺序获取基准点相邻的点云为待测点，判断该待测点是否为断裂点，若为断裂点则保留该点云；若为非断裂点，则以该待测点构建有向向量并计算向量夹角ρ，若ρ不大于阈值，则保留该待测点否则滤除该待测点；具体公式为：

式中，B为基准点,P_i为待测点；O为扫描中心点；

(4.2)若上述待测点被保留，则将该保留的待测点替换为基准点，重复步骤(4.1)直至遍历扫描线上所有的点云；若上述待测点被滤除，则不替换该基准点B，重复步骤(4.1)直至遍历扫描线上所有的点云；

(5)采用二次曲面拟合法滤除测量区域内属于非地面点的点云，获取最终的滤波后的点云。

2.根据权利要求1所述基于断裂线约束的机载点云快速滤波方法，其特征在于，步骤(1)中所述提取断裂点具体包括以下步骤：

(1.1)计算测量区域中每个点云与其邻域内点云的平均曲率值H；具体公式为：

式中，E＝r_x·r_x、F＝r_x·r_y、G＝r_y·r_y、

r_x，r_y，r_xx，r_xy，r_yy均为z(x,y)的偏导数；z(x,y)＝ax²+bxy+cy²为任意一点云及其邻域内的点云使用最小二乘法拟合曲面并对其进行参数化后的表达式，a、b、c均为拟合参数；

(1.2)针对测量区域中每个点云，计算高程小于该点的邻域点数与全部邻域点数的比值G；

(1.3)根据各点云的H和G计算各点云的特征值D，保留特征值D大于预设值的点云作为断裂点；具体公式为：

D＝|H|q₁+Gq₂

式中，q₁为Ⅰ类断裂点权重，q₂为Ⅱ类断裂点权重。

3.根据权利要求1所述基于断裂线约束的机载点云快速滤波方法，其特征在于，步骤(3)中所述的粗过滤具体包括以下步骤：

(3.1)定义每条扫描线上最左端的点云为初始点P₁，按照扫描线上的点云顺序选取与初始点相邻的点云，计算该点云与初始点的高度差，若高度差小于设定高度差阈值，则该点云为地面点记为P₂；否则滤除该点继续按照扫描线上的点云顺序选取后续新点重复上述操作，直至确定地面点P₂；

(3.2)选取与P_i+1相邻的点云记为P_i+2，以P_i P_i+1、P_i+1P_i+2长度分别构建等边三角形P_i O_iP_i+1、等边三角形P_i+1O_i+1P_i+2；若∠O_iP_i+1O_i+1小于夹角阈值，则判定P_i+2点为地面点；若∠O_iP_{i+ 1}O_i+1大于夹角阈值，则判定该P_i+2点为非地面点，滤除该P_i+2点并继续按照扫描线上的点云顺序选择后续新点作为新的P_i+2，重新构建等边三角形并获取夹角∠O_iP_i+1O_i+1，再将∠O_iP_{i+ 1}O_i+1与夹角阈值比较直至获取为地面点的P_i+2点；i取1……n，n为正整数；

(3.3)令i自加1，重复步骤(3.2)直至后续新点为扫描线上最后一个点云，滤除扫描线上的非地面点，确定扫描线上所有地面点；

(3.4)计算扫描线上每个点云的点云平坦度，保留点云平坦度大于平坦度阈值的点云。

4.根据权利要求1所述基于断裂线约束的机载点云快速滤波方法，其特征在于，步骤(5)具体包括以下步骤：

(5.1)根据原始点云数据中每个点云的横纵坐标范围划分格网，不同网格之间形成重叠区域、非重叠区域；获取落入单元格内的滤波后扫描线上的m个地面点，采用最小二乘原理求得误差平方和σ²达到最小时的单元格的曲面参数A′、B′、C′、D′、E′、F′具体的公式为：

式中，z(x_i,y_i)为点云的坐标；

(5.2)确定每个格网曲面的参数后，判断原始点云中每个点云所在的格网位置：若原始点云落入格网非重叠区域，则计算该点云在单个格网的拟合高程z(x,y)；若原始点云落入格网重叠区域，则计算该点云在不同格网曲面拟合高程z(x,y)的平均值；z(x,y)获取的公式为：

z(x,y)＝A'x²+B'y²+C'xy+D'x+E'y+F'

(5.3)比较各点云的拟合高程值与实际高程的差值，滤除差值大于预设阈值的点云，获得最终的滤波后的点云。

5.根据权利要求3所述基于断裂线约束的机载点云快速滤波方法，其特征在于，步骤(3)中所述的粗过滤还包括：

(3.5)对于点云中的任一点云，计算并统计高程大于该点的邻域点数与全部邻域点数的比值，保留比值小于设定阈值的点云。

6.根据权利要求3所述基于断裂线约束的机载点云快速滤波方法，其特征在于，步骤(3.2)中由等边三角形几何关系计算获得O_i的坐标定义具体公式为：

式中，P_i表示为

P_i+1表示为

P_i+2表示为(x_i+2 ^p，y_i+2 ^p，z_i+2 ^p)，O₁表示为(x_i ^o，y_i ^o，z_i ^o)。

7.根据权利要求3所述基于断裂线约束的机载点云快速滤波方法，其特征在于，步骤(3.4)中点云平坦度的计算公式为：

式中，X_i表示为点云的领域，0＜i≤n；D表示为协方差矩阵；

为空间邻域的中心点；λ为点云平坦度；λ₁、λ₂、λ₃分别为协方差矩阵D在三维空间的特征值。

8.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至权利要求7任一所述方法的步骤。

9.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至权利要求7任一所述方法的步骤。

Images