CN111598780A - 一种适用于机载LiDAR点云的地形自适应插值滤波方法 - Google Patents
一种适用于机载LiDAR点云的地形自适应插值滤波方法 Download PDFInfo
- Publication number
- CN111598780A CN111598780A CN202010407153.3A CN202010407153A CN111598780A CN 111598780 A CN111598780 A CN 111598780A CN 202010407153 A CN202010407153 A CN 202010407153A CN 111598780 A CN111598780 A CN 111598780A
- Authority
- CN
- China
- Prior art keywords
- point
- points
- ground
- grid
- elevation
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 132
- 238000001914 filtration Methods 0.000 title claims abstract description 97
- 230000003044 adaptive effect Effects 0.000 title claims abstract description 19
- 230000002159 abnormal effect Effects 0.000 claims abstract description 33
- 230000008030 elimination Effects 0.000 claims abstract description 6
- 238000003379 elimination reaction Methods 0.000 claims abstract description 6
- 239000013256 coordination polymer Substances 0.000 claims description 24
- 230000008569 process Effects 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 14
- 238000011160 research Methods 0.000 claims description 7
- 238000012795 verification Methods 0.000 claims description 3
- 230000000750 progressive effect Effects 0.000 abstract description 7
- 230000009286 beneficial effect Effects 0.000 abstract description 3
- 238000007689 inspection Methods 0.000 abstract 1
- 238000011439 discrete element method Methods 0.000 description 22
- 230000006870 function Effects 0.000 description 11
- 230000008901 benefit Effects 0.000 description 8
- 238000012545 processing Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 6
- 230000000877 morphologic effect Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 230000000903 blocking effect Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000010801 machine learning Methods 0.000 description 4
- ATJFFYVFTNAWJD-UHFFFAOYSA-N Tin Chemical compound [Sn] ATJFFYVFTNAWJD-UHFFFAOYSA-N 0.000 description 2
- 238000013459 approach Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000010187 selection method Methods 0.000 description 1
- 238000012706 support-vector machine Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4023—Scaling of whole images or parts thereof, e.g. expanding or contracting based on decimating pixels or lines of pixels; based on inserting pixels or lines of pixels
-
- 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
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/194—Segmentation; Edge detection involving foreground-background segmentation
-
- 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/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10044—Radar image
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种适用于机载LiDAR点云的地形自适应插值滤波方法。该方法包括低异常点去除、地面种子点选取、地面种子点核查以及地面点渐进选择等步骤。其中，在低异常点去除步骤中，本发明采用高程直方图实现低异常点自动剔除，以降低人工干预程度；在地面种子点选取步骤中，本发明利用一维离散光滑样条获取大量地面种子点，以提高初始地面参考面精度；在地面点渐进选择步骤中，本发明采用尺度无关插值估计待分类点值，以避免空间位置误差影响，此外，本发明还采用地形自适应滤波阈值以适应不同地形场景点云，以降低不同地形场景区点云错分和漏分的水平。本发明方法利于提升点云滤波的精度。
Description
技术领域
本发明涉及一种适用于机载LiDAR点云的地形自适应插值滤波方法。
背景技术
目前，机载激光雷达(LiDAR)技术以其主动性和高效性等优势已成为三维空间数据采集的主要手段，广泛应用于数字高程模型生产、三维城市构建、森林资源调查以及山体滑坡监测等领域。区别于直接接触式数据采集方式，原始机载LiDAR点云数据包括各种地面点和非地面点，因此，原始机载LiDAR点云数据在使用前必须进行滤波处理。
过去二十年间，国内外众多研究人员提出了各种点云滤波算法。根据工作原理，现有的滤波算法可分为坡度滤波、形态学滤波、插值滤波、分块滤波以及机器学习滤波等。
坡度滤波假设相邻地面点坡度较小，因此，如果两点间坡度超过设定阈值，则较高的点为非地面点。然而，该滤波方法对坡度阈值非常敏感，特别在地形复杂区域容易将地形特征点错分为非地面点。为此，相关人员提出了一系列自适应坡度滤波方法，确保坡度阈值随地形复杂度自适应变化。由于坡度滤波仅利用相邻点之间的局部信息，该方法在地形断裂区和陡坡区滤波误差较大。
形态学滤波主要借助形态学开运算对点云滤波，即以一定尺寸的结构元开运算栅格化点云，如果处理前后栅格点高程差超过规定阈值，则对应点标记为非地面点。其中，结构元尺寸显著影响滤波精度：该尺寸必须足够大以能够去除各种尺寸非地面物，同时，该尺寸必须足够小以尽可能保留地形细节信息。为此，相关人员提出了一系列渐进形态学滤波算法以剔除各种尺寸非地面物。然而，形态学滤波以栅格化点云为操作对象，因此在点云栅格化过程中不可避免的造成精度损失。
插值滤波首先对地面种子点插值构建初始地面参考面，接着比较待滤波点与地面参考面的距离，如果该距离小于设定阈值，则对应的点标记为地面点，最后用更新后的地面种子点重复上述过程直至没有点标记为地面点。
作为插值滤波的典型代表，渐进三角网滤波因被纳入商业软件TerraSolid而被广泛使用。然而，该滤波方法使用TIN插值构建地面参考面，因此在地形复杂区域滤波精度较低。
为此，部分学者建议使用较高精度的非线性插值方法(如线性估计、薄板样条和克里金等)代替TIN插值，并提出了多分辨分层插值滤波方法。该类滤波方法在计算过程中需要构建一系列不同分辨率DEM作为地面参考面，并比较每个待分类点到对应DEM网格点的距离。然而，当待分类点与DEM网格中心点不重合时，势必产生空间位置误差，而且该空间位置误差与地形坡度呈正比，即坡度越大，误差越大。因此，插值滤波在陡坡区域滤波精度较低。
分块滤波以对点云分割后的块为基本处理单元，并利用块和块之间的拓扑关系、块和地面参考面的高程差等对各个块滤波。鉴于点云块可捕捉到地形断裂点，理论上分块滤波在地形断裂区具有较好的滤波效果。然而，分块滤波的优势是以对点云准确分块为前提，不准确分块容易抵消分块滤波的高效性。而且，该类滤波在地形复杂区域由于各个块中点数较少，优势并不明显。
近年来，机器学习方法被广泛使用点云滤波。常用的机器学习算法包括人工神经网络、条件随机场、支持向量机和深度学习等。然而，机器学习滤波方法需要大量的训练样本，制作过程费时费力；而且，该类滤波具有计算量庞大以及人工参与程度高等缺陷。
根据上述分析，表1总结了上述五种滤波方法的优缺点。
表1各种机载LiDAR点云滤波方法优缺点比较
大量实验分析表明，现有的各种滤波方法在场景简单区可以取得较好的滤波精度，但在场景复杂区域容易引起点云漏分或错分问题，降低了滤波精度，后续需要大量人工编辑工作，费时费力，而且点云编辑精度完全取决于操作人员的熟练程度，主观性强。
发明内容
本发明的目的在于提出一种适用于机载LiDAR点云的地形自适应插值滤波方法，以提升各种场景下点云滤波的精度。本发明为了实现上述目的，采用如下技术方案：
一种适用于机载LiDAR点云的地形自适应插值滤波方法，包括如下步骤：
I.利用高程直方图自动去除原始点云中的低异常点；
II.选取地面种子点，具体过程如下：
首先将原始点云的研究区域覆盖一定尺寸的网格，并选择每个网格中的最低点作为候选点；记录每个候选点的三维坐标，而且候选点按照栅格方式存储；
然后以栅格为操作对象，分别逐行和逐列选择初始地面种子点；
最后将行方式和列方式都选中的点标记为地面种子点；
其中，从栅格的任意一行或列栅格数据中选择初始地面种子点的步骤为：
II.1.初始地面种子点选择；
利用一维移动窗口选择初始地面种子点；其中，选中的网格点含有高程值，未选中的网格点标记为NaN；一维移动窗口的尺寸取值为最大建筑物的短边长，设为S；
II.2.样条插值；
借助一维离散光滑样条法对初始地面种子点插值，估计其它未选中网格点的值；
II.3.网格点高程差计算；
首先利用原始网格点高程减去插值后网格点高程，得到网格点高程差；然后将网格点高程差小于设定阈值eth的点标记为地面点；
II.4.地面种子点更新；
利用步骤II.3中选中的地面点更新地面种子点；
II.5.重复上述步骤II.2-步骤II.4，直至该行或列栅格数据没有地面点被选中；
III.利用稳健曲面拟合方法对步骤II中选取的地面种子点进行核查，将非地面点剔除；
IV.对地面点进行渐进选择；
经过以上步骤I-步骤III中的低异常点剔除、地面种子点选择以及核查后，原始点云被分为三部分，即非地面点集NonGP、地面点集GP以及候选点集CP；
其中，非地面点集NonGP中仅包括低异常点，地面点集GP中包括地面种子点，候选点集CP中既包括地面点也包括非地面点；
下面从候选点集CP中选择地面点，候选点集CP中剩余的点即为非地面点集NonGP；
以薄板样条TPS作为插值函数，对地面点渐进选择的具体步骤为：
IV.1.对于候选点集中任一待分类点P，利用kd树从地面点集GP选择点P的k个邻近点；
IV.2.以步骤IV.1中k个邻近点为已知点，利用TPS插值求算点P高程值；
IV.3.对点P的观测值和估算值做差，求算点P的插值误差e；
IV.4.如果插值误差e小于阈值t，则点P标记为地面点；
阈值t的计算公式为：t＝t0+c*slope；
其中，t0为初始误差阈值，c为坡度系数，slope为点P处地形坡度值；
利用薄板样条TPS基本方程求算slope的值，计算公式为：
其中：
式中，fx表示点P处函数值在x方向偏导数，fy表示点P处函数值在y方向偏导数；
αj表示第j个基函数系数；
(x,y)表示点P的平面坐标，(xj,yj)表示点P的第j个邻近地面种子点平面坐标；
rj表示点P与第j个邻近地面种子点的欧式距离，ε为大于0的常数；
IV.5.重复上述步骤IV.1-步骤IV.4，直至候选点集CP中所有待分类点分类完毕；
IV.6.将选中的地面点从候选点集CP中删除，利用选中的地面点更新地面点集GP；
IV.7.重复上述步骤IV.1-步骤IV.6，直至候选点集CP中没有地面点被选中，则候选点集CP中剩余的点标记为非地面点集NonGP。
优选地，步骤I具体为：
I.1.将原始点云的高程值分段，即将整个高程值的范围从左到右均匀分成m个间隔Ci，i＝1,2,…,m，每个间隔的高程跨度L＝(Hmax-Hmin)/m；
其中，Hmax表示原始点云的高程最大值，Hmin表示原始点云的高程最小值；
I.2.计算每个间隔Ci中包含的点数；
I.3.如果第一个间隔C1中点数超过设定阈值n，则原始点云中认定没有低异常点；否则，从间隔C1到间隔Ck中选择点数小于n的最右间隔，并将该最右间隔记为Cr；
I.4.计算最右间隔Cr中包含的点中高程最大值MAX，并将该高程该最大值MAX作为区分原始点云中正常点和低异常点的阈值；
将高程值大于MAX的点自动分类为正常点，否则自动分类为低异常点。
优选地，步骤II中，在原始点云研究区域覆盖的网格尺寸等于原始点云中点平均距离；
其中，点平均距离是指所有相邻点距离的平均值。
优选地，步骤III中，利用稳健曲面拟合方法对地面种子点进行核查的具体过程如下：
对于步骤II选取的地面种子点中任一点R，利用kd树选中点R周围的k个邻近点，并对这些邻近点进行平面拟合，然后计算点R到该拟合平面的距离；
如果该距离大于设定阈值rth，则该点R为非地面点，应从地面种子点中剔除。
优选地，步骤III中，在进行平面拟合时，利用最小截断二乘法拟合切平面；
最小截断二乘法的目标函数为：
式中，(r2)(i)表示第i个邻近点到拟合平面的距离平方，且(r2)(1)≤(r2)(2)≤…≤(r2)(k)，h＝k/2。
本发明具有如下优点：
如上所述，本发明提出了一种用于机载LiDAR点云的地形自适应插值滤波方法，该方法利用高程直方图自动去除低异常点，以降低人工交互工作量；借助一维离散样条插值方法选取尽量多的地面种子点，提升初始地面参考面精度；直接利用已选取的地面种子点估算待分类点的高程，避免使用各种尺度DEM带来的空间位置误差影响；构建自适应滤波阈值，以降低不同地形场景区点云错分和漏分的水平。本发明利于提升不同场景下点云滤波的精度。
附图说明
图1为本发明实施例中适用于机载LiDAR点云的地形自适应插值滤波方法的流程框图。
图2为本发明实施例中基于高程直方图的低异常点剔除示意图。
图3为本发明实施例中按照行方式选择地面种子点的流程示意图。
图4为本发明实施例中最小截断二乘法与传统最小二乘法拟合平面的结果比较示意图。
图5为本发明实施例中对地面点渐进选择的流程示意图。
图6为空间位置误差对滤波结果影响示意图。
图7为高密度点云区正射影像图。
图8为高密度点云DEM山体阴影图。
图9为本发明方法和PMF、PTD、MHF滤波方法在样区3的DEM示意图。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
如图1所示，一种适用于机载LiDAR点云的地形自适应插值滤波方法，包括如下步骤：
I.利用高程直方图自动去除原始点云中的低异常点。
将原始点云中的低异常点剔除，以免错选为初始地面种子点，受LiDAR系统误差以及多路径影响，原始点云中不可避免的含有异常值。异常值分两类：低异常点和高异常点。
前者明显低于其它点云，而后者明显高于其它点云。
传统的插值滤波方法在计算过程中，需要选择局部邻域最低点作为初始地面种子点，因此低异常点容易被错选。本实施例采用高程直方图自动去除原始点云中的低异常点。
以上低异常点自动去除方式具有客观性和高效性，有效降低了低异常点被错选的概率。
如图2所示，利用高程直方图自动去除原始点云中的低异常点的具体过程为：
I.1.将原始点云的高程值分段，即将整个高程值的范围从左到右均匀分成m个间隔Ci，i＝1,2,…,m，每个间隔的高程跨度L＝(Hmax-Hmin)/m。
其中，Hmax表示原始点云的高程最大值，Hmin表示原始点云的高程最小值。
I.2.计算每个间隔Ci中包含的点数。
I.3.如果第一个间隔C1中点数超过设定阈值n，则原始点云中认定没有低异常点；否则，从间隔C1到间隔Ck中选择点数小于n的最右间隔，并将该最右间隔记为Cr。
I.4.计算最右间隔Cr中包含的点中高程最大值MAX，并将该高程该最大值MAX作为区分原始点云中正常点和低异常点的阈值。
将高程值大于MAX的点自动分类为正常点，否则自动分类为低异常点。
在本实施例中m例如取10，n例如取值为50，k＝m/2。
II.选取地面种子点。
对插值滤波而言，地面种子点的数量和质量显著影响滤波精度。
而传统插值滤波方法主要借助局部最小值法选择地面种子点，即使用稍大于研究区最大建筑物尺寸的移动窗口按一定步长逐步移动，并选择每个窗口中的最低点作为地面种子点。
该方法虽然能保证种子点质量，但是获取的点数量较少，难以保证初始地面参考面精度。
基于此，本实施例构建了一种一维离散光滑样条选点法，以保证地面种子点数量和质量。该方法主要从栅格化点云中选择地面种子点，即
首先将原始点云的研究区域覆盖一定尺寸的网格，选择每个网格中的最低点作为候选点；
在原始点云研究区域覆盖的网格尺寸等于原始点云中点平均距离；
其中，该点平均距离是指所有相邻点距离的平均值。
候选点按照两种方式保存：(1)逐行记录每个候选点的三维坐标，(2)按照候选点所在的行列值以栅格方式存储，标记为G。
然后以栅格G为操作对象，分别逐行和逐列选择初始地面种子点，如图1所示。
最后将行方式和列方式都选中的点标记为地面种子点。
如图3所示，以行方式选点为例，从栅格G中任意一行的栅格数据(如图3(a)所示的栅格候选点(Z))中选择地面种子点步骤为：
II.1.初始地面种子点选择，如图3(b)所示。
本实施例利用一维移动窗口选择初始地面种子点。其中，选中的网格点含有高程值，未选中的网格点标记为NaN。一维移动窗口的尺寸取值为最大建筑物的短边长，设为S。
II.2.样条插值，如图3(c)所示。
本实施例借助一维离散光滑样条法对初始地面种子点插值，估计其它未选中网格点的值。
该一维离散光滑样条法具有计算速度快、低内存需求以及自动化程度高等优势，因此能够保证地面种子点选择的准确性和高效性。
II.3.网格点高程差计算，如图3(d)所示。
首先利用原始网格点高程Z减去插值后网格点高程Zs，得到网格点高程差(Z-Zs)；然后将网格点高程差(Z-Zs)小于设定阈值eth的点标记为地面点。
本实施例中设定阈值eth例如取值为0.5m。
II.4.地面种子点更新，如图3(e)所示。
本实施例利用步骤II.3中选中的地面点更新地面种子点。
II.5.重复上述步骤II.2-步骤II.4，直至该行栅格数据没有地面点被选中。
至此完成了从栅格G中任意一行的栅格数据中选择地面种子点。对于从栅格G中其他行的栅格数据以及任意一列的栅格数据中选择地面种子点的过程与上面过程相同。
需要说明的是，在图3中出现的ID移动窗口，表示1维移动窗口。
III.按照上述方式选择的地面种子点不可避免的含有非地面点。为了提升地面种子点的质量，本实施例借助稳健的曲面拟合方法对选取的地面种子点核查，将非地面点剔除。
利用曲面拟合方法对选取的地面种子点进行核查的具体过程如下：
对于步骤II选取的地面种子点中任一点R，利用kd树选中点R周围的k(例如k＝12)个邻近点，并对这些邻近点进行平面拟合，然后计算点R到该拟合平面的距离；
若距离大于设定阈值rth(例如rth＝0.1m)，则该点R为非地面点，从地面种子点中剔除。
为了提升曲面拟合的稳健性，本实施例利用最小截断二乘法(Least TrimmedSquares method,LTS)拟合切平面；最小截断二乘法LTS的目标函数为：
式中，(r2)(i)表示第i个邻近点到拟合平面的距离平方，且(r2)(1)≤(r2)(2)≤…≤(r2)(k)，h＝k/2。
图4给出了利用最小截断二乘法与传统最小二乘法拟合平面的结果比较示意图。
由图4可见，如果用经典的最小二乘(Least Squares method,LS)拟合切平面，第1、2、8和9点均被错分为非地面点；而以LTS拟合切平面，第7点被准确标记为非地面点。
因此，相比于LS拟合方式，LTS拟合方式能够准确捕捉到地面种子点中的非地面点。
IV.对地面点进行渐进选择。
经过以上步骤I-步骤III中的低异常点剔除、地面种子点选择以及核查后，原始点云被分为三部分，即非地面点集NonGP、地面点集GP以及候选点集CP。
其中，非地面点集NonGP中仅包括低异常点，地面点集GP中包括地面种子点，候选点集CP中既包括地面点也包括非地面点。
下面从候选点集CP中选择地面点，候选点集CP中剩余的点即为非地面点集NonGP。
本实施例采用的插值方法为薄板样条(Thin Plate Spline,TPS)，其基本方程为：
式中，f(x,y)为待分类点P的高程值，αj表示第j个基函数系数；
(x,y)表示待分类点P的平面坐标，k为待分类点P周围邻近的地面种子点数量；
rj表示待分类点P与第j个邻近地面种子点的欧式距离。
q(x)为径向基函数，对薄板样条TPS而言，q(x)＝x2ln(x2)；
α和b表示薄板样条TPS的未知系数，它们需要通过求算如下方程获得，即：
式中，(Q)ij＝q(rij)；e＝(1 … 1)T；α＝(α1 … αk)T；f＝(f1 … fk)T；
rij表示第i个点和第j个已知点的欧式距离；(Q)ij表示矩阵Q的第i行j列元素值；
fj表示第j个邻近点地面点的高程值，j的取值范围为[1,k]。
如图5所示，以TPS作为插值函数，则本实施例对地面点渐进选择的步骤为：
IV.1.对于候选点集中任一待分类点P，利用kd树从地面点集GP选择点P的k个邻近点。
IV.2.以步骤IV.1中k个邻近点为已知点，利用TPS插值求算点P高程值，具体计算公式如公式(2)和公式(3)所示。本实施例通过直接利用已选取的地面种子点估算待分类点P的高程，能够避免使用各种尺度DEM带来的空间位置误差影响。
IV.3.对点P的观测值和估算值做差，求算点P的插值误差e。
IV.4.如果插值误差e小于阈值t，则点P标记为地面点。
鉴于TPS插值在地形复杂区域插值误差要大于平坦区域，因此，阈值t应随着地形复杂度变化，如图1所示，即在地形平坦区域较小，而在地形复杂区域较大。
本实施例中阈值t的计算公式为：t＝t0+c*slope；
其中，t0为初始误差阈值，c为坡度系数，slope为点P处地形坡度值。
为了进一步避免空间位置误差影响，本实施例利用薄板样条TPS基本方程(如上面公式(2))求算slope的值，计算公式为：
其中：
式中，fx表示点P处函数值在x方向偏导数，fy表示点P处函数值在y方向偏导数；
(xj,yj)表示点P的第j个邻近地面种子点平面坐标；
ε为防止计算公式中分母为0而人为加的较小的常数，例如ε＝0.001。
IV.5.重复上述步骤IV.1-步骤IV.4，直至候选点集CP中所有待分类点分类完毕。
IV.6.将选中的地面点从候选点集CP中删除，利用选中的地面点更新地面点集GP。
IV.7.重复上述步骤IV.1-步骤IV.6，直至候选点集CP中没有地面点被选中，则候选点集CP中剩余的点标记为非地面点集NonGP。
区别于以往借助不同分辨率DEM估计待分类点值，本实施例直接利用TPS对待分类点插值，进而避免了DEM网格中心点和待求点位置不重合而引起的空间位置误差影响。
例如，由图6可见，在网格A和D处，地形平坦，空间位置误差对滤波结果影响较小，所有点准确分类，而在网格B和C处，地形坡度较大，受空间位置误差影响，部分点出现错分。因此，尺度无关空间插值能够有效避免空间位置误差影响，利于提升后续滤波精度。
此外，为了适应地形复杂度变化的不同，本实施例还提出了地形自适应滤波阈值t，即在地形平坦区域，坡度slope较小，使得t值较小；而在陡坡区域，坡度slope较大，使得t值偏大，由此能够准确捕捉到地形特征点，进而降低滤波误差影响。
下面结合具体实例对本发明方法的有效性进行说明：
首先以ISPRS提供的15组基准数据为数据源，分析本发明方法的计算精度，然后将滤波结果与近10年来(2010-2020)提出的滤波方法进行比较。
表2给出ISPRS基准数据的基本统计信息，主要包括地面点数(#BE)，非地面点数(#OBJ)，点密度、平均坡度以及主要场景特征。
表2 ISPRS15组数据的基本统计信息
鉴于ISPRS数据于二十年前发布，受到当时数据采集设备软硬件影响，点密度较低(0.16-1.02pts/m2)。为验证本发明方法处理高密度点云的可行性，本实施例进一步以四组城区和林区高密度机载LiDAR点云为数据源，分析其滤波精度。
此外，本实施例还将本发明计算结果与现有的三种经典滤波方法比较，包括渐进形态学滤波(PMF)、渐进三角网加密滤波(PTD)以及多尺度层次插值滤波(MHF)。
表3列出了四组点云数据的基本统计信息。
除此之外，图7和8分别给出了它们的正射影像图和DEM山体阴影图。
其中，图7(a)-图7(d)分别表示样区1至样区4的高密度点云区正射影像图，图8(a)-图8(d)分别表示样区1至样区4的高密度点云区DEM山体阴影图。
由此可见，林区(样区3和4)地形复杂度明显高于城区(样区1和2)。
表3四组高密度机载LiDAR点云基本统计信息
本发明实施例主要借助I类误差(T.I.)、II类误差(T.II.)、总误差(T.E)以及kappa系数(κ)分析本发明方法和其它滤波方法的精度。
本发明包括三个主要参数：一维移动窗口尺寸(S)，初始高差阈值(t0)和坡度系数(c)。类似于其它滤波方法，本发明针对不同样区的滤波结果均是在参数最优条件下获得。
表4列出了本发明方法对15组基准数据的最优参数以及对应的滤波误差。
整体而言，本发明方法在15组样区中的平均总误差仅为2.70％，kappa系数高达90.84％，表明本发明提出的新滤波方法在不同场景区具有较好的滤波效果。
就单个样区而言，本发明方法在样区31精度最高，总误差和kappa系数分别为0.78％和98.44％。这是因为该样区地形平坦、场景简单，如表2所示。
由于样区11建筑物和低矮植被主要分布在陡坡上，本发明方法在该区总误差最大，为8.72％；本发明方法在样区53kappa系数最小，为68.35％，这主要归因于该区地面点数明显少于非地面点数，使得少量错分地面点导致较大I类误差，进而拉低kappa系数。
表4本发明方法对ISPRS基准数据最优参数以及滤波误差
表5给出了本发明方法与三种经典滤波方法对高密度点云数据处理结果。
就I类误差而言，PMF在样区1、2、3精度高于其它方法，而本发明方法在样区4精度最优；以II类误差为精度指标，本发明方法在样区1和3结果精度最高，在样区2和4仅次于最优方法PTD和MHF。以总误差和kappa系数为精度指标，本发明方法在样区1、3、4精度明显优于其它方法，且在样区2精度与最优方法PMF接近。
总体而言，本发明方法的平均总误差均小于其它方法，而kappa系数均高于其它方法。
具体而言，本发明方法总误差比PMF、PTD、MHF分别小11.4％、66.4％、23.3％；而kappa系数比PMF、PTD、MHF分别高1.6％、23.8％and 2.9％。
表5本发明方法与经典滤波方法对高密度点云数据处理结果比较(％)
鉴于利用滤波后地面点构建的DEM精度可间接反映各种方法的滤波效果，本实施例计算了本发明方法和其它三种经典滤波方法在各个样区DEM中误差。
表6表明所有方法在城区的DEM精度优于林区，这是因为前者的地形复杂度要高于后者，如图8所示。除了样区2外，本发明方法中误差明显小于其它三种方法，以平均中误差为精度指标，本发明方法精度分别比PMF、PTD、MHF高3.7％、50.4％、48.1％。
由此证明，本发明方法处理高密度机载LiDAR点云的有效性。
表6本发明方法与经典方法在高密度点云区DEM中误差(m)
图9展示了本发明方法和其它三种滤波方法在样区3的DEM示意图。结果表明，各种方法构建的DEM均不同程度受滤波误差影响，其中，PTD结果最差，DEM表面非常粗糙，表明大量低矮植被错分为地面点，如图3d所示。PMF构建的DEM尽管明显优于PTD，但在地形细节处不够光滑，少量非植被点被错分为地面点，如图3b所示。MHF在左上角处也有少些低矮植被错分为地面点，如图3c所示。整体而言，尽管本发明方法也受到I类误差影响，图3a所示，但它构建的DEM与参考DEM最接近，精度也最高。
由以上分析可见，无论是处理低密度的基准数据还是高密度的机载LiDAR点云数据，本发明方法均表现了一定的优势，这主要归功于本发明方法中所采用的自动低异常点剔除、大量地面种子点获取、尺度无关插值以及地形自适应滤波阈值等技术手段。
表7展示了低异常点自动剔除技术手段监测到的低异常数、该异常点中包含的非地面点数以及在有无使用该技术手段条件下本发明方法滤波精度。
结果表明，除了样区12外，该技术手段均能准确捕捉到样区中的低异常点。在有无使用该技术手段条件下，除了样区41外，本发明方法总误差并没有明显变化。这是因为：
一方面，各个样区低异常点较少，它们对滤波精度影响甚微；另一方面，本发明方法在地面种子点选择过程中，借助一维光滑样条插值以及采用LTS拟合切平面，这两种手段均具有较高的抗差性，由此可有效抑制低异常点影响，使得不被错选为地面种子点。
表7低异常点自动剔除技术手段精度分析
表8进一步展示了其它技术手段在不同样区对本发明方法精度影响程度。
结果表明，对地面种子点获取而言，以样区12为例，使用该技术手段后本发明方法精度均有一定程度的上升，其总误差比没有使用该技术手段降低8.2％，kappa系数提升0.49％。
表8不同技术手段对本发明方法滤波精度影响
就尺度无关插值而言，以样区23为例，使用该技术手段后本发明方法精度上升明显，其总误差比没有使用该技术手段降低17.2％，kappa系数提升1.9％。由实验表明，在没有使用尺度无关插值条件下，大量地形断裂点被错分为非地面点，导致I类误差明显上升。这是因为，相同尺度DEM表达不同复杂度地形时，平坦区精度较高，而在复杂区精度较低。
就自适应滤波阈值而言，以样区52为例，使用该技术手段后本发明精度有大幅度升高，其总误差比没有使用该技术手段降低54.5％，kappa系数提升13.9％。由实验表明，在没有使用自适应滤波阈值条件下，陡坡区大量地面点被错分为非地面点，导致I类误差非常明显。
通过以上实例分析，进一步验证了本发明方法所采用的各种技术手段的有效性。综上，本发明方法所采用的低异常点剔除、地面种子点加密、尺度无关插值和坡度自适应滤波阈值等技术手段均在一定程度上提升了滤波精度，使其有效适应不同场景下点云滤波。
当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。
Claims (5)
1.一种适用于机载LiDAR点云的地形自适应插值滤波方法，其特征在于，
包括如下步骤：
I.利用高程直方图自动去除原始点云中的低异常点；
II.选取地面种子点，具体过程如下：
首先将原始点云的研究区域覆盖一定尺寸的网格，并选择每个网格中的最低点作为候选点；记录每个候选点的三维坐标，而且候选点按照栅格方式存储；
然后以栅格为操作对象，分别逐行和逐列选择初始地面种子点；
最后将行方式和列方式都选中的点标记为地面种子点；
其中，从栅格的任意一行或列栅格数据中选择初始地面种子点的步骤为：
II.1.初始地面种子点选择；
利用一维移动窗口选择初始地面种子点；其中，选中的网格点含有高程值，未选中的网格点标记为NaN；一维移动窗口的尺寸取值为最大建筑物的短边长，设为S；
II.2.样条插值；
借助一维离散光滑样条法对初始地面种子点插值，估计其它未选中网格点的值；
II.3.网格点高程差计算；
首先利用原始网格点高程减去插值后网格点高程，得到网格点高程差；然后将网格点高程差小于设定阈值eth的点标记为地面点；
II.4.地面种子点更新；
利用步骤II.3中选中的地面点更新地面种子点；
II.5.重复上述步骤II.2-步骤II.4，直至该行或列栅格数据没有地面点被选中；
III.利用稳健曲面拟合方法对步骤II中选取的地面种子点进行核查，将非地面点剔除；
IV.对地面点进行渐进选择；
经过以上步骤I-步骤III中的低异常点剔除、地面种子点选择以及核查后，原始点云被分为三部分，即非地面点集NonGP、地面点集GP以及候选点集CP；
其中，非地面点集NonGP中仅包括低异常点，地面点集GP中包括地面种子点，候选点集CP中既包括地面点也包括非地面点；
下面从候选点集CP中选择地面点，候选点集CP中剩余的点即为非地面点集NonGP；
以薄板样条TPS作为插值函数，对地面点渐进选择的具体步骤为：
IV.1.对于候选点集中任一待分类点P，利用kd树从地面点集GP选择点P的k个邻近点；
IV.2.以步骤IV.1中k个邻近点为已知点，利用TPS插值求算点P高程值；
IV.3.对点P的观测值和估算值做差，求算点P的插值误差e；
IV.4.如果插值误差e小于阈值t，则点P标记为地面点；
阈值t的计算公式为：t＝t0+c*slope；
其中，t0为初始误差阈值，c为坡度系数，slope为点P处地形坡度值；
利用薄板样条TPS基本方程求算slope的值，计算公式为：
其中：
式中，fx表示点P处函数值在x方向偏导数，fy表示点P处函数值在y方向偏导数；
αj表示第j个基函数系数；
(x,y)表示点P的平面坐标，(xj,yj)表示点P的第j个邻近地面种子点平面坐标；
rj表示点P与第j个邻近地面种子点的欧式距离，ε为大于0的常数；
IV.5.重复上述步骤IV.1-步骤IV.4，直至候选点集CP中所有待分类点分类完毕；
IV.6.将选中的地面点从候选点集CP中删除，利用选中的地面点更新地面点集GP；
IV.7.重复上述步骤IV.1-步骤IV.6，直至候选点集CP中没有地面点被选中，则候选点集CP中剩余的点标记为非地面点集NonGP。
2.根据权利要求1所述的地形自适应插值滤波方法，其特征在于，
所述步骤I具体为：
I.1.将原始点云的高程值分段，即将整个高程值的范围从左到右均匀分成m个间隔Ci，i＝1,2,…,m，每个间隔的高程跨度L＝(Hmax-Hmin)/m；
其中，Hmax表示原始点云的高程最大值，Hmin表示原始点云的高程最小值；
I.2.计算每个间隔Ci中包含的点数；
I.3.如果第一个间隔C1中点数超过设定阈值n，则原始点云中认定没有低异常点；否则，从间隔C1到间隔Ck中选择点数小于n的最右间隔，并将该最右间隔记为Cr；
I.4.计算最右间隔Cr中包含的点中高程最大值MAX，并将该高程该最大值MAX作为区分原始点云中正常点和低异常点的阈值；
将高程值大于MAX的点自动分类为正常点，否则自动分类为低异常点。
3.根据权利要求1所述的地形自适应插值滤波方法，其特征在于，
所述步骤II中，在原始点云研究区域覆盖的网格尺寸等于原始点云中点平均距离；
其中，点平均距离是指所有相邻两点间距离的平均值。
4.根据权利要求1所述的地形自适应插值滤波方法，其特征在于，
所述步骤III中，利用稳健曲面拟合方法对地面种子点进行核查的具体过程如下：
对于步骤II选取的地面种子点中任一点R，利用kd树选中点R周围的k个邻近点，并对这些邻近点进行平面拟合，然后计算点R到该拟合平面的距离；
如果该距离大于设定阈值rth，则该点R为非地面点，应从地面种子点中剔除。
5.根据权利要求4所述的地形自适应插值滤波方法，其特征在于，
所述步骤III中，在进行平面拟合时，利用最小截断二乘法拟合切平面；
最小截断二乘法的目标函数为：
式中，(r2)(i)表示第i个邻近点到拟合平面的距离平方，且(r2)(1)≤(r2)(2)≤…≤(r2)(k)，h＝k/2。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010407153.3A CN111598780B (zh) | 2020-05-14 | 2020-05-14 | 一种适用于机载LiDAR点云的地形自适应插值滤波方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010407153.3A CN111598780B (zh) | 2020-05-14 | 2020-05-14 | 一种适用于机载LiDAR点云的地形自适应插值滤波方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111598780A true CN111598780A (zh) | 2020-08-28 |
| CN111598780B CN111598780B (zh) | 2022-03-18 |
Family
ID=72187391
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010407153.3A Active CN111598780B (zh) | 2020-05-14 | 2020-05-14 | 一种适用于机载LiDAR点云的地形自适应插值滤波方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111598780B (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113589319A (zh) * | 2021-10-08 | 2021-11-02 | 中南大学 | 一种迭代极小值激光雷达点云滤波方法 |
| CN114612806A (zh) * | 2022-02-07 | 2022-06-10 | 安徽理工大学 | 一种提高消费级无人机dem产品精度方法 |
| CN116579949A (zh) * | 2023-05-31 | 2023-08-11 | 浙江省测绘科学技术研究院 | 适用于城市多噪声环境下机载点云地面点滤波方法 |
| CN117194928A (zh) * | 2023-11-07 | 2023-12-08 | 湖南中云图地理信息科技有限公司 | 基于gnss的地理形变监测系统 |
| WO2024060209A1 (zh) * | 2022-09-23 | 2024-03-28 | 深圳市速腾聚创科技有限公司 | 一种处理点云的方法和雷达 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106897686A (zh) * | 2017-02-19 | 2017-06-27 | 北京林业大学 | 一种机载lidar电力巡检点云分类方法 |
| US20180081035A1 (en) * | 2016-09-22 | 2018-03-22 | Beijing Greenvalley Technology Co., Ltd. | Method and device for filtering point cloud data |
| CN110335352A (zh) * | 2019-07-04 | 2019-10-15 | 山东科技大学 | 一种机载激光雷达点云的双基元多分辨率层次滤波方法 |
-
2020
- 2020-05-14 CN CN202010407153.3A patent/CN111598780B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20180081035A1 (en) * | 2016-09-22 | 2018-03-22 | Beijing Greenvalley Technology Co., Ltd. | Method and device for filtering point cloud data |
| CN106897686A (zh) * | 2017-02-19 | 2017-06-27 | 北京林业大学 | 一种机载lidar电力巡检点云分类方法 |
| CN110335352A (zh) * | 2019-07-04 | 2019-10-15 | 山东科技大学 | 一种机载激光雷达点云的双基元多分辨率层次滤波方法 |
Non-Patent Citations (2)
| Title |
|---|
| CHUANFA CHEN,MENGYING WANG,BINGTAO CHANG,YANYAN LI: "Multi-Level Interpolation-Based Filter", 《IEEE ACCESS》 * |
| 李乐林等: "一种顾及地形复杂度的LiDAR点云多尺度滤波方法", 《测绘科学》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113589319A (zh) * | 2021-10-08 | 2021-11-02 | 中南大学 | 一种迭代极小值激光雷达点云滤波方法 |
| CN114612806A (zh) * | 2022-02-07 | 2022-06-10 | 安徽理工大学 | 一种提高消费级无人机dem产品精度方法 |
| WO2024060209A1 (zh) * | 2022-09-23 | 2024-03-28 | 深圳市速腾聚创科技有限公司 | 一种处理点云的方法和雷达 |
| CN116579949A (zh) * | 2023-05-31 | 2023-08-11 | 浙江省测绘科学技术研究院 | 适用于城市多噪声环境下机载点云地面点滤波方法 |
| CN117194928A (zh) * | 2023-11-07 | 2023-12-08 | 湖南中云图地理信息科技有限公司 | 基于gnss的地理形变监测系统 |
| CN117194928B (zh) * | 2023-11-07 | 2024-01-26 | 湖南中云图地理信息科技有限公司 | 基于gnss的地理形变监测系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111598780B (zh) | 2022-03-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111598780B (zh) | 一种适用于机载LiDAR点云的地形自适应插值滤波方法 | |
| CN106529469B (zh) | 基于自适应坡度的无人机载LiDAR点云滤波方法 | |
| CN110992341A (zh) | 一种基于分割的机载LiDAR点云建筑物提取方法 | |
| CN103077529A (zh) | 基于图像扫描的植物叶片特征分析系统 | |
| CN111340723B (zh) | 一种地形自适应的机载LiDAR点云正则化薄板样条插值滤波方法 | |
| CN110335352B (zh) | 一种机载激光雷达点云的双基元多分辨率层次滤波方法 | |
| CN111667574B (zh) | 从倾斜摄影模型自动重建建筑物规则立面三维模型的方法 | |
| CN114266987B (zh) | 一种无人机高边坡危岩体智能识别方法 | |
| CN112669333A (zh) | 一种单木信息提取方法 | |
| CN116523898A (zh) | 一种基于三维点云的烟草表型性状提取方法 | |
| CN113177897A (zh) | 一种无序3d点云的快速无损滤波方法 | |
| CN112308966A (zh) | 基于多级曲率约束的高斯混合模型点云滤波方法 | |
| CN114862715A (zh) | 一种融合地形特征语义信息的tin渐进加密去噪方法 | |
| CN109242786B (zh) | 一种适用于城市区域的自动化形态学滤波方法 | |
| CN116186864B (zh) | 一种基于bim技术的深基坑模型快速建模方法及系统 | |
| CN114743059B (zh) | 一种综合地形地貌特征的海底地理实体自动分类方法 | |
| CN112381029B (zh) | 一种基于欧氏距离的机载LiDAR数据建筑物提取方法 | |
| CN115410036A (zh) | 一种高压架空输电线路关键要素激光点云自动分类方法 | |
| CN114283163A (zh) | 一种多平台激光雷达大田玉米高通量茎叶分离方法 | |
| CN112215184A (zh) | 一种基于三维激光扫描仪的油茶果树产量检测方法 | |
| CN117132478B (zh) | 一种基于法向量二范数特征参数的轨道点云去噪方法 | |
| Zou et al. | An adaptive strips method for extraction buildings from light detection and ranging data | |
| CN116704333B (zh) | 一种基于激光点云数据的单株林木检测方法 | |
| CN115294485B (zh) | 一种市政工程测量定位方法及系统 | |
| CN117152172A (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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |