一种地形自适应的规定密度机载LiDAR点云精简方法

摘要

本发明涉及一种地形自适应的规定密度机载LiDAR点云精简方法。首先对原始点云划分虚拟格网，并根据格网内地形起伏度和粗糙度，将格网划分成平坦区、粗糙区和非粗糙区三种不同的地表类型，然后对不同地表类型区域实行不同的抽稀策略，得到规定密度的点云。本发明针对不同地表类型区域采用不同的抽稀策略，能充分保留地形起伏关键点，防止局部密度过低和空洞的产生，避免出现关键点缺失导致的地面模型的失真。本发明参数设置简单，参数数量较少，且具有一定的地形自适应能力，可有效解决实际生产中同一区域内不同地表形态点云数据混杂导致的阈值设置困难的问题。实验证明本方法处理大数据量的点云数据速度快、精度高，能满足实际生产需求。

权利要求

1.一种地形自适应的规定密度机载LiDAR点云精简方法，其特征在于，包括如下步骤：步骤1，建立二维格网点云索引；

步骤2，计算格网内地形起伏度；

步骤3，根据地形起伏度阈值α将格网划分为平坦区和非平坦区；

步骤4，针对非平坦区，计算格网内地形粗糙度；

地形粗糙度的计算方法如下：

式中，Rough为地形粗糙度，Zi为第i个点的实测高程，n为局部范围内的点个数，Z(xi,yi)为第i个点在二次曲面上的拟合高程，二次曲面可由格网内局部点云数据进行二次拟合得到，计算方式如下：

2 2

Z＝a0+a1x+a2y+a3x+a4y+a5xy (3)式中，Z表示局部区域内点的高程，x和y为局部区域内点的平面坐标，a0、a1、a2、a3、a4和a5为二次曲面函数的系数；

对于数据冗余的格网区域，采用最小二乘法求解，即拟合高程与实测高程误差的平方和最小，公式如下：式中，Z(xi,yi)为第i个点的拟合高程，Zi为第i个点的实测高程，n为局部范围内的点个数；

建立误差方程，并求解二次曲面函数的系数矩阵，公式如下：V＝BX‑L (5)

其中，

T ‑1 T

X＝(BB) BL (6)

式中，V表示拟合高程与实测高程之间的误差矩阵，B表示公式(3)中各系数对应的坐标组合项，X为公式(3)中待求的二次曲面函数的系数矩阵，L为实测高程矩阵；

步骤5，根据粗糙度阈值β，将非平坦地区进一步划分为粗糙区和非粗糙区；

步骤6，针对平坦区地表格网，依据指定密度计算出格网所需点数，采取随机降采样获得抽稀后的点云数据；

步骤7，针对非平坦区点云，计算格网内点云的曲率，保留高程极值点和高曲率点；

步骤8，根据步骤7保留的点数与规定密度下单位格网所需点数量的大小关系，对非平坦区不同区域采用不同的抽稀策略。

2.如权利要求1所述的一种地形自适应的规定密度机载LiDAR点云精简方法，其特征在于：所述步骤1是根据点云数据的XY坐标将点云投影到二维平面，并划分规则格网，用格网的行列号实现点云数据的组织与管理。

3.如权利要求1所述的一种地形自适应的规定密度机载LiDAR点云精简方法，其特征在于：所述步骤2中地形起伏度是指局部范围的高程极大值与极小值之差，反映了局部范围的地形起伏，地形起伏度的计算方法如下：dH＝Hmax‑Hmin (1)

式中，dH为局部范围的地形起伏度，Hmax、Hmin分别为局部范围内地面点的高程极大值和极小值。

4.如权利要求1所述的一种地形自适应的规定密度机载LiDAR点云精简方法，其特征在于：所述步骤3中是将格网内地形起伏度大于α的区域标记为平坦区，地形起伏度小于α的标记为非平坦区，α为经验值。

5.如权利要求1所述的一种地形自适应的规定密度机载LiDAR点云精简方法，其特征在于：所述步骤5中是将格网内地形粗糙度大于β的区域标记为粗糙区，地形粗糙度小于β的区域标记为非粗糙区，β为经验值。

6.如权利要求1所述的一种地形自适应的规定密度机载LiDAR点云精简方法，其特征在于：所述步骤7中局部高程极值点包括格网内高程极大值点和极小值点两个，可以利用根据离散点高程直接算出；高曲率点是曲率值大于阈值γ的点，γ为经验值，假设目标点Pi＝T(xi,yi,zi)的协方差矩阵为CPi，则：

式中， 为目标点及其邻域点的中心位置，N为目标点的邻域点个数，CPi为对称半正定矩阵；

曲率curve的计算方式如下：

式中，λ1、λ2、λ3为协方差矩阵的三个特征值，且λ1≥λ2≥λ3。

7.如权利要求1所述的一种地形自适应的规定密度机载LiDAR点云精简方法，其特征在于：所述步骤8中若高程极值点和高曲率点保留的数量多于规定密度下单位格网所需点数量，则对非粗糙区，按照曲率从大到小的顺序保留至规定的数量高曲率点和高程极值点，对粗糙区保留全部曲率点和高程极值点；若高程极值点和高曲率点数量少于规定密度下单位格网所需点数量，则计算与所需点数之间的点数差，在剩余点中以随机采样方式补齐，处理完全部点云数据后获得抽稀过后的点云。

说明书

一种地形自适应的规定密度机载LiDAR点云精简方法

技术领域

[0001] 本发明属于LiDAR数据处理与应用领域，特别是涉及一种地形自适应的规定密度机载LiDAR点云精简方法。

背景技术

[0002] 随着硬件技术的不断进步，获取的LiDAR点云数据密度已可达到每平方米几十个点。高密度的点云虽然有助于更细化地描述地形特征，但对于较平坦区域却带来了大量的数据冗余，这些冗余数据将直接影响数据存储、处理和显示速度的效率，使现有点云处理算法、软件受到挑战，同时也提高了对硬件的要求。点云抽稀是在确保满足生产精度或实际需求的基础上，对原始点云进行有规则地筛选，最大限度地精简点云数量。根据不同算法采用的规则，点云抽稀可归纳为随机采样算法和顾及地形特征的算法。随机采样算法主要包括基于规则格网的抽稀算法与基于系统的抽稀算法。前者将整个数据范围划分为很多的小格网，随机选择格网中一个数据点或几个数据点的平均值代替格网中的所有数据。后者是在加载原始点云时，设置采样间距N，取N个点云数据的第一个或随机从这N个点云中抽取一个数据点。两个算法同样简单快捷，但均不能很好地保留地形特征，降低了点云数据的精度。顾及地形特征的算法中，多采用高程、高差、曲率、坡度等参数。常见的有Pamelas提出的DDR(Data Density Reduction)抽稀算法、徐景中等人提出的基于点云离散度的抽稀算法、嫪志修等人提出的基于坡度的抽稀算法、杨明军等人提出的约束TIN节点抽稀法和侯文广等人提出在测地空间中采用泊松碟采样的抽稀算法等。在基于点云数据的测绘产品生产过程中，如何在减小点云数据量降低冗余的同时保留点云数据中的关键点、使点云数据精度达标成为研究的重点。

发明内容

[0003] 本发明针对现有技术的不足，提供一种地形自适应的规定密度机载LiDAR点云精简方法。首先对原始点云划分虚拟格网，计算格网内地形起伏度，并依据地形起伏度阈值将格网划分为平坦区和非平坦区，然后对非平坦区格网计算格网地形粗糙度，根据设置的阈值进一步将其划分为粗糙区和非粗糙区，针对平坦区，采取指定密度的随机抽稀方法；针对非平坦区，计算格网内点云曲率，保留高曲率点和高程极值点，并根据保留的高曲率点和高程极值点数量与规定密度下单位格网所需点数量的大小关系，对非平坦区不同区域采用不同的抽稀策略，得到规定密度的点云。

[0004] 为了达到上述目的，本发明提供的技术方案是一种地形自适应的规定密度机载LiDAR点云精简方法，包括以下步骤：

[0005] 步骤1，建立二维格网点云索引；

[0006] 步骤2，计算格网内地形起伏度；

[0007] 步骤3，根据地形起伏度阈值α将格网划分为平坦区和非平坦区；

[0008] 步骤4，针对非平坦区，计算格网内地形粗糙度；

[0009] 步骤5，根据粗糙度阈值β，将非平坦地区进一步划分为粗糙区和非粗糙区；

[0010] 步骤6，针对平坦区地表格网，依据指定密度计算出格网所需点数，采取随机降采样获得抽稀后的点云数据；

[0011] 步骤7，针对非平坦区点云，计算格网内点云的曲率，保留高程极值点和高曲率点；

[0012] 步骤8，根据步骤7保留的点数与规定密度下单位格网所需点数量的大小关系，对非平坦区不同区域采用不同的抽稀策略。

[0013] 而且，所述步骤1中是根据点云数据的XY坐标将点云投影到二维平面，并划分规则格网，用格网的行列号实现点云数据的组织与管理。

[0014] 而且，所述步骤2中地形起伏度是指局部范围的高程极大值与极小值之差，反映了局部范围的地形起伏，地形起伏度的计算方法如下：

[0015] dH＝Hmax‑Hmin (1)

[0016] 式中，dH为局部范围的地形起伏度，Hmax、Hmin分别为局部范围内地面点的高程极大值和极小值。

[0017] 而且，所述步骤3中是将格网内地形起伏度大于α的区域标记为平坦区，地形起伏度小于α的标记为非平坦区，α为经验值。

[0018] 而且，所述步骤4中地形粗糙度表达的是局部地形的粗糙程度或者褶皱程度，以局部范围内各地面点的高程拟合误差为基础计算地形粗糙度，计算方法如下：

[0019]

[0020] 式中，Rough为地形粗糙度，Zi为第i个点的实测高程，n为局部范围内的点个数，Z(xi,yi)为第i个点在二次曲面上的拟合高程，二次曲面可由格网内局部点云数据进行二次拟合得到，计算方式如下：

[0021] Z＝a0+a1x+a2y+a3x2+a4y2+a5xy (3)

[0022] 式中，Z表示局部区域内点的高程，x和y为局部区域内点的平面坐标，a0、a1、a2、a3、a4和a5为二次曲面函数的系数。

[0023] 对于数据冗余的格网区域，包含的点数一般大于6个，存在多余观测，因此采用最小二乘法求解，即拟合高程与实测高程误差的平方和最小，公式如下：

[0024]

[0025] 式中，Z(xi,yi)为第i个点的拟合高程，Zi为第i个点的实测高程，n为局部范围内的点个数。

[0026] 建立误差方程，并求解二次曲面函数的系数矩阵，公式如下：

[0027] V＝BX‑L (5)

[0028] 其中，

[0029] X＝(BTB)‑1BTL (6)

[0030] 式中，V表示拟合高程与实测高程之间的误差矩阵，B表示公式(3)中各系数对应的坐标组合项，X为公式(3)中待求的二次曲面函数的系数矩阵，L为实测高程矩阵。

[0031] 而且，所述步骤5中是将格网内地形粗糙度大于β的区域标记为粗糙区，地形粗糙度小于β的区域标记为非粗糙区，β为经验值。

[0032] 而且，所述步骤7中局部高程极值点包括格网内高程极大值点和极小值点两个，可以利用根据离散点高程直接算出，高曲率点是曲率值大于阈值γ的点，γ为经验值，假设目T标点Pi＝(xi,yi,zi)的协方差矩阵为CPi，则：

[0033]

[0034]

[0035] 式中， 为目标点及其邻域点的中心位置，N为目标点的邻域点个数，CPi为对称半正定矩阵。

[0036] 对协方差矩阵进行特征值分解，得到三个特征值λ1、λ2、λ3(λ1≥λ2≥λ3)，及其对应的特征向量e1、e2、e3。e3为估算的法向量，λ3描述了曲面沿法向量方向的变化量，曲率curve的计算方式如下：

[0037]

[0038] 虽然基于特征值估算的曲率的精度略低，但是处理速度更快，对于大范围冗余的机载LiDAR地面点云数据具有更好的适用性，因此采用公式(9)估算激光脚点的曲率特征，邻域类型选择K邻域。

[0039] 而且，所述步骤8中若高程极值点和高曲率点保留的数量多于规定密度下单位格网所需点数量，则对非粗糙区，按照曲率从大到小的顺序保留至规定的数量高曲率点和高程极值点，对粗糙区保留全部曲率点和高程极值点。若高程极值点和高曲率点数量少于规定密度下单位格网所需点数量，则计算与所需点数之间的点数差，在剩余点中以随机采样方式补齐。处理完全部点云数据后获得抽稀过后的点云。

[0040] 与现有技术相比，本发明具有如下优点：

[0041] (1)准确性：本方法针对不同地形区域采用不同的抽稀策略，充分保留地形起伏关键点，防止局部密度过低和空洞的产生，同时也避免了关键点缺失导致的地面模型的失真；

[0042] (2)易用性：参数设置简单，参数数量较少，且具有一定的地形自适应能力，可以有效解决实际生产中同一片区域内不同地表形态点云数据混杂导致的阈值设置困难的问题；

[0043] (3)高效性，实验证明本方法针对大数据量的点云数据处理速度快，精度高，可以满足实际生产需求。

附图说明

[0044] 图1是本发明实施例的技术流程图。

[0045] 图2是本发明实施例山地地区数据不同方法的实验结果对比图，其中图2(a)是原始点云数据，图2(b)是系统抽稀方法抽稀后的点云数据，图2(c)是曲率特征全局方法抽稀后的点云数据，图2(d)是本发明方法抽稀后的点云数据。

[0046] 图3是本发明实施例山地地区数据的实验结果局部细节对比图，其中图3(a)是系统抽稀方法抽稀后的点云数据局部细节图，图3(b)是本发明方法抽稀后的点云数据局部细节图。

[0047] 图4是本发明实施例在城市地区数据1的实验结果局部细节对比图，其中图4(a)是原始点云数据，图4(b)是本发明方法抽稀后的点云数据。

[0048] 图5是本发明实施例城市地区数据数据2的实验结果局部细节对比图，其中图5(a)是原始点云数据，图5(b)是本发明方法抽稀后的点云数据。

[0049] 图6是本发明实施例丘陵地区数据1的实验结果局部细节对比图，其中图6(a)是原始点云数据，图6(b)是本发明方法抽稀后的点云数据。

[0050] 图7是本发明实施例丘陵地区数据1的实验结果局部细节对比图，其中图7(a)是原始点云数据，图7(b)是本发明方法抽稀后的点云数据。

[0051] 图8是本发明实施例丘陵地区数据2的实验结果局部细节对比图，其中图8(a)是原始点云数据，图8(b)是本发明方法抽稀后的点云数据。

[0052] 图9是本发明实施例丘陵地区数据2的实验结果局部细节对比图，其中图9(a)是原始点云数据，图9(b)是本发明方法抽稀后的点云数据。

[0053] 图10是本发明实施例丘陵地区数据3的实验结果局部细节对比图，其中图10(a)是原始点云数据，图10(b)是本发明方法抽稀后的点云数据。

具体实施方式

[0054] 本发明针对现有技术的不足，提供一种地形自适应的规定密度机载LiDAR点云精简方法。首先对原始点云划分虚拟格网，计算格网内地形起伏度，并依据地形起伏度阈值将格网划分为平坦区和非平坦区，然后对非平坦区格网计算格网地形粗糙度，根据设置的阈值进一步将其划分为粗糙区和非粗糙区，针对平坦区，采取指定密度的随机抽稀方法；针对非平坦区，计算格网内点云曲率，保留高曲率点和高程极值点，并根据保留的高曲率点和高程极值点数量与规定密度下单位格网所需点数量的大小关系，对非平坦区不同区域采用不同的抽稀策略，得到规定密度的点云。

[0055] 下面结合附图和实施例对本发明的技术方案作进一步说明。

[0056] 如图1所示，本发明实施例的流程包括以下步骤：

[0057] 步骤1，建立二维格网点云索引。根据点云数据的XY坐标将点云投影到二维平面，并划分规则格网，用格网的行列号实现点云数据的组织与管理。

[0058] 步骤2，计算格网内地形起伏度。地形起伏度是指局部范围的高程极大值与极小值之差，反映了局部范围的地形起伏，地形起伏度的计算方法如下：

[0059] dH＝Hmax‑Hmin (1)

[0060] 式中，dH为局部范围的地形起伏度，Hmax、Hmin分别为局部范围内地面点的高程极大值和极小值。

[0061] 步骤3，根据地形起伏度阈值α将格网划分为平坦区和非平坦区。若格网内地形起伏度大于α，则将其标记为平坦区，否则标记为非平坦区，α为经验值。

[0062] 步骤4，针对非平坦区，计算格网地形粗糙度。地形粗糙度表达了局部地形的粗糙程度或者褶皱程度，以局部范围内各地面点的高程拟合误差为基础计算地形粗糙度，计算方法如下：

[0063]

[0064] 式中，Rough为地形粗糙度，Zi为第i个点的实测高程，n为局部范围内的点个数，Z(xi,yi)为第i个点在二次曲面上的拟合高程，二次曲面可由格网内局部点云数据进行二次拟合得到，计算方式如下：

[0065] Z＝a0+a1x+a2y+a3x2+a4y2+a5xy (3)

[0066] 式中，Z表示局部区域内点的高程，x和y为局部区域内点的平面坐标，a0、a1、a2、a3、a4和a5为二次曲面函数的系数。

[0067] 对于数据冗余的格网区域，包含的点数一般大于6个，存在多余观测，因此采用最小二乘法求解，即拟合高程与实测高程误差的平方和最小，公式如下：

[0068]

[0069] 式中，Z(xi,yi)为第i个点的拟合高程，Zi为第i个点的实测高程，n为局部范围内的点个数。

[0070] 建立误差方程，并求解二次曲面函数的系数矩阵，公式如下：

[0071] V＝BX‑L (5)

[0072] 其中，

[0073] X＝(BTB)‑1BTL (6)

[0074] 式中，V表示拟合高程与实测高程之间的误差矩阵，B表示公式(3)中各系数对应的坐标组合项，X为公式(3)中待求的二次曲面函数的系数矩阵，L为实测高程矩阵。

[0075] 步骤5，根据粗糙度阈值β将非平坦地区进一步划分为粗糙区和非粗糙区。若格网内地形粗糙度大于β，则将其标记为粗糙区，否则标记为非粗糙区，β为经验值。

[0076] 步骤6，针对平坦区地表格网，依据指定密度计算出格网所需点数，采取随机降采样获得抽稀后的点云数据。

[0077] 步骤7，针对非平坦区点云，计算格网内部点云的曲率，保留高程极值点和高曲率点。

[0078] 局部高程极值点包括格网内高程极大值点和极小值点两个，反映了局部区域高程的变化范围，可以利用根据离散点高程直接算出。

[0079] 点的曲率大小反映了地形曲面在该点处的弯曲程度，因此高曲率点是表达局部地形发生明显变化的重要特征点，高曲率点的集合能够刻画整个区域的关键地形骨架。将曲T率值大于阈值γ的点作为高曲率点，假设目标点Pi＝(xi,yi,zi)的协方差矩阵为CPi，则：

[0080]

[0081]

[0082] 式中， 为目标点及其邻域点的中心位置，N为目标点的邻域点个数，CPi为对称半正定矩阵。

[0083] 对协方差矩阵进行特征值分解，得到三个特征值λ1、λ2、λ3(λ1≥λ2≥λ3)，及其对应的特征向量分别为e1、e2、e3。e3为估算的法向量，λ3描述了曲面沿法向量方向的变化量，曲率curve计算方式如下：

[0084]

[0085] 虽然基于特征值估算的曲率的精度略低，但是处理速度更快，对于大范围冗余的机载LiDAR地面点云数据具有更好的适用性。因此采用公式(9)估算激光脚点的曲率特征，邻域类型选择K邻域。

[0086] 步骤8，根据步骤7保留的点数与规定密度下单位格网所需点数量的大小关系，对非平坦区不同区域采用不同的抽稀策略。

[0087] 若高程极值点和高曲率点保留的数量多于规定密度下单位格网所需点数量，则对非粗糙区，按照曲率从大到小的顺序保留至规定的数量高曲率点和高程极值点，对粗糙区保留全部曲率点和高程极值点。若高程极值点和高曲率点数量少于规定密度下单位格网所需点数量，则计算与所需点数之间的点数差，在剩余点中以随机采样方式补齐。处理完全部点云数据后获得抽稀过后的点云。

[0088] 图2是不同抽稀方案的对比结果图，系统抽稀方法的结果虽然基本保留了总体地形特征，但部分山谷区域的特征表现并不明显；基于曲率特征进行全局抽稀的方法虽然保留了很多地形特征点，但出现了大范围的数据空洞；相比之下，本发明提出的抽稀方法不仅较好地保留了地形特征，特别是部分山谷区域的特征经过抽稀后更加明显，而且避免了数据空洞的出现，得到的抽稀效果最好。

[0089] 图3是系统抽稀方法和本发明方法抽稀后的点云数据局部细节对比图，通过观察局部细节可以发现，本发明的抽稀方法将山谷区域的点识别为重要地形特征点，并进行了针对性的保留，而系统抽稀方法是间隔性采样，对局部地形特征点的保留比较随机和盲目。在相同的点云保留率的条件下，本发明的抽稀方法更能兼顾地形特征，保留关键地形特征点的能力更强，点位分布也更为合理。

[0090] 图4‑10为本发明在各类点云数据的实验结果，可以看到均取得了良好的实验效果，不论是城区的建筑物轮廓点还是山地的山脊、山谷，丘陵等地形特征区域关键点都得到了很好地保留，同时点云数量得到大幅度的降低，极大地减少了数据的冗余，证明了本发明的泛用性和实用性。

[0091] 本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。