CN103729846B - 基于不规则三角网的LiDAR点云数据边缘检测方法 - Google Patents

Abstract

本发明公开了一种基于不规则三角网的LiDAR点云数据边缘检测方法：读取目标区域原始LiDAR点云数据；剔除粗差噪声点；将剔除粗差噪声后的LiDAR点云数据进行Delaunay三角剖分；计算不规则三角网中每个空间三角形的形变量并保存；计算三角形的形变量阈值；将得到的形变量大于形变量阈值的空间三角形标记为边缘三角形；计算狭长三角形；比较各三角形的狭长度L_i和狭长度阈值；确定边缘点，得到LiDAR点云数据的边缘图像。本发明的方法利用三角剖分过程中三角形的形状变化进行具有高程突变的边缘点检测，利用狭长三角形进行数据空白区域的边缘点检测，实现了在原始数据中进行边缘检测保留LiDAR点云数据更多的信息。

Description

基于不规则三角网的LiDAR点云数据边缘检测方法

技术领域

本发明属于遥感应用研究领域，具体涉及一种基于不规则三角网的LiDAR点云数据边缘检测方法。

背景技术

机载LiDAR***能够直接获取地面三维数据，具有高精度、高密度、高效率和低成本的优点，利用它快速获得的海量数据可得到所需的各种图像产品，因此在现代测绘中发挥越来越重要的作用。图像边缘是图像的重要特性之一，是图像中比较有价值的信息。将检测到的边缘信息作为特征信息，为图像的配准与融合提供了前提条件。因此，利用机载LiDAR***快速获取高精度的LiDAR点云数据，并检测出理想的边缘具有重要的意义。LiDAR点云数据是指机载LiDAR***通过发射和接收激光脉冲得到地表的高精度3维点坐标。

传统的LiDAR点云数据的边缘检测方法，一般是将离散的LiDAR点云数据进行规则格网重采样，得到数字深度(高程)矩阵，再将不同的深度值进行灰度量化和拉伸变换，得到同灰度图像一样的距离灰度图像，这个过程称为“距离成像”。赖旭东、万幼川等对深度图像进行中值滤波，利用Robert算子、Sobel算子以及Prewitt算子进行边缘提取；王大莹、程新文、潘慧波等对深度图像使用最佳阈值进行二值化后使用Canny和Log算子进行边缘检测；徐景中、万幼川等对深度图像利用自适应的中值滤波以及自适应阈值的Canny算子进行边缘检测；吴杭彬、刘春等对深度图像利用数学形态学中的膨胀和腐蚀方法进行序贯运算。并将得到的图像进行边缘提取和边缘矢量化,得到每个地物所对应的边缘和数据点。

通过对上述各种LiDAR点云数据的边缘检测方法进行分析，发现它们存在有如下缺陷：距离成像过程中进行重采样时会带来误差，使得一些河流或者点云空白区域的边缘无法被检测出来；在灰度量化和拉伸变换时会带来舍入误差，高程差异较小的点在进行灰度量化时可能会量化为同一个灰度级或者灰度差异很小，也可能将差别较小的高程量化为不同的灰度级，因此在边缘检测时都会带来一定的误差，使得一些边缘点检测不到或得到错误的边缘点。

综上，现有的基于深度图像进行LiDAR点云数据的边缘检测方法虽然能够达到一定的效果，但是深度图像化的过程中却丢失了一些数据信息，得到的边缘信息不够完整，因此得到的边缘信息并不能得到LiDAR点云数据的其他等有用信息，不利于后续的处理。

发明内容

针对现有的将LiDAR点云数据转换为深度图像进行边缘检测的不足，本文的目的在于，提供一种基于不规则三角网的LiDAR点云数据边缘检测方法，该方法直接在原始数据上基于不规则三角网中三角形形变量的差异情况进行LiDAR点云数据边缘检测。该方法充分利用三角剖分过程中三角形的形状变化进行具有高程突变的边缘点检测，利用狭长三角形进行数据空白区域的边缘点检测，实现了在原始数据中进行边缘检测保留LiDAR点云数据更多的信息。本发明可以应用于机载LiADR点云数据边缘检测。

为了实现上述任务，本发明采用如下技术方案予以解决：

一种基于不规则三角网的LiDAR点云数据边缘检测方法，具体包括如下步骤：

步骤1：读取目标区域原始LiDAR点云数据；

步骤2：剔除粗差噪声点；

步骤3：将步骤2得到的剔除粗差噪声后的LiDAR点云数据进行Delaunay三角剖分，生成xoy平面上的不规则三角网；在不规则三角网上加上各点的高程信息后形成xyz空间不规则三角网TIN；

步骤4：计算不规则三角网中每个空间三角形的形变量并保存；

步骤5：计算三角形的形变量阈值；

步骤6：将步骤4得到的形变量大于形变量阈值的空间三角形标记为边缘三角形；如果形变量大于形变量阈值,进入步骤8，否则进入步骤7；

步骤7：计算狭长三角形，比较各三角形的狭长度L_i和狭长度阈值Thresd，如果L_i大于Thresd，则该三角形为狭长三角形，将该三角形标记为边缘三角形，其中，Thresd不小于100；

步骤8：确定边缘点，得到LiDAR点云数据的边缘图像，具体步骤如下：

801：粗边缘点的确定：

对于步骤6得到的每个边缘三角形利用式8计算高程阈值H(G)，并分别对每个边缘三角形中的顶点进行粗边缘点的选取：

$H\left(G\right)=\frac{H\left(A\right)+H\left(B\right)+H\left(C\right)}{3}---\left(8\right)$

其中，H(A)、H(B)、H(C)分别是A、B、C点的高程信息，单位：m；H(G)为高程阈值，单位：m；

若边缘三角形的三个顶点A,B,C中有任意两点A,B的高程均大于H(G)且|H(A)-H(B)|≤σ，σ取0.5-2.5，则将A,B两点均取为粗边缘点，否则，将A,B,C三点中高程最大的点取为粗边缘点；

对步骤7得到的每个边缘三角形的三个顶点都确定为粗边缘点；

:802：删除所有粗边缘点中的孤立点，将剩余的粗边缘点作为最终的边缘点输出，即得到LiDAR点云数据的边缘图像。

进一步的，所述步骤2的粗差噪声点利用式1计算得到：

n＝∑(|P_i-P_j|＜σ) (1)

其中，P_i，P_j表示粗边缘点中点i和点j的坐标，n表示原始LiDAR点云数据中与坐标P_i的距离小于σ的点的个数，n取0-10，σ取10-20m。

进一步的，所述步骤4的计算不规则三角网中每个空间三角形的形变量是指：对每个三角形利用式2进行计算：

$s=C\×\underset{i=1}{\overset{3}{\Σ}}|\sin A_i^{\′}-\sin A_i|---\left(2\right)$

其中，s是空间三角形的形变量，C是常数，取大于1的任意值，A_i表示xoy平面三角形的内角，A′_i表示添加高程信息形成空间三角形后与A_i对应的内角，i＝1,2,3。

进一步的，所述步骤5计算三角形的形变量阈值的具体步骤如下：

501：读取所有空间三角形的形变量s组成的形变量集合S；设定形变量阈值T^k的初始值T⁰：

$T^0=\frac{S_{\min }+S_{max}}{2}---\left(3\right)$

其中，S_min和S_max分别表示形变量集合S中的最小值和最大值；

502：将形变量集合S中所有的三角形形变量s分成形变量大于T^k的集合A和形变量小于等于T^k的集合B两部分；

503：利用式4、式5计算A、B两部分的平均形变量g_A和g_B：

$g_A=\frac{1}{n}\underset{S\left(i\right)>T^k}{\Σ}S\left(i\right)---\left(4\right)$

$g_B=\frac{1}{m}\underset{S\left(i\right)\≤T^k}{\Σ}S\left(i\right)---\left(5\right)$

其中，S(i)是形变量集合S中的第i个三角形形变量，n和m分别是A和B中三角形形变量s的个数；

504：利用式6计算新的阈值T^k+1：

$T^{k+1}=\frac{g_A+g_B}{2}---\left(6\right)$

505：判断条件|T^k+1-T^k|≤δ,k＝0,1,2,…是否成立，δ是允许误差，此处误差范围C*0.01-C*0.05，其中，C是式2中的常数；条件成立，则形变量阈值T^k更新为步骤504得到的T^k+1，然后执行步骤6；否则，阈值仍为T^k，此时令k＝k+1，并返回步骤503。

进一步的，所述步骤7的计算狭长三角形是指：

利用式7计算步骤3得到的空间三角网中所有三角形的狭长度；

L_i＝|A_iG_i|²+|B_iG_i|²+|C_iG_i|² (7)

其中，L_i表示三角网中第i个三角形的狭长度；A_i，B_i，C_i分别为第i个三角形的三个顶点，G_i为第i个三角形的重心。

进一步的，所述步骤8的802中所述孤立点利用步骤2的公式9进行计算：

n＝∑(|P_i-P_j|＜σ) (9)

其中，P_i表示粗边缘点中点i的XY坐标，n表示粗边缘点中与坐标P_i的距离小于σ的点的个数，n取0-10，σ取10-20m。

与现有技术相比，本发明具有以下优点：

第一，本发明在对具有高程突变的地物边缘时，应用Delaunay三角剖分建立不规则三角网，利用三角形的形变量大小作为边缘三角形的判断准则，将形变量大于阈值的三角形标记位边缘三角形，对边缘三角形处理得到边缘点。该方法无需将LiDAR点云数据生成深度图像，从而避免现有的边缘检测方法中在转化为深度图像时的边缘数据损失，提高了检测到边缘的精度，也有利于后续的LiDAR数据特征提取以及与遥感图像和SAR图像的融合。

第二，本发明在对有河流等数据空白区域等LiDAR点云数据进行边缘检测时，利用狭长三角形作为测度，对空白区域边缘进行检测，消除了数据格网化过程中数据丢失以及内插造成的边缘信息丢失，保持了边缘信息的完整性。

附图说明

图1为本发明的基于不规则三角网的LiDAR点云数据边缘检测方法的总流程图。

图2为本发明的方法用于LiDAR点云数据边缘检测的仿真试验1的效果图。其中，图2(a)是原始LiDAR点云数据图像；图2(b)是对原始图像进行三角剖分后的图像；图2(c)是本发明的方法得到的边缘点检测图；图2(d)是将原始数据和边缘点叠加的效果图。

图3为本发明的方法用于含有河流等数据空白区域的数据边缘检测的仿真试验2的效果图。其中，图3(a)是原始LiDAR点云数据图像；图3(b)是对原始图像进行三角剖分后的图像；图3(c)是本发明的方法得到的边缘点检测图；图3(d)是将原始数据和边缘点叠加的效果图。

以下结合附图和具体实施方式对本发明进一步解释说明。

具体实施方式

参照图1，本发明的基于不规则三角网的LiDAR点云数据边缘检测方法，具体包括如下步骤：

步骤1：读取目标区域原始LiDAR点云数据；

步骤2：剔除粗差噪声；

原始点云数据会产生一些粗差噪声，为了得到准确的点云数据，首先要对粗差噪声点进行剔除，粗差噪声点利用式1确定：

n＝∑(|P_i-P_j|＜σ) (1)

其中，P_i，P_j表示粗边缘点中点i和点j的坐标，n表示原始LiDAR点云数据中与坐标P_i的距离小于σ的点的个数，n取0-10，σ取10-20m；

步骤3：将步骤2得到的剔除粗差噪声后的LiDAR点云数据进行Delaunay三角剖分，即将不规则分布的离散数据点生成xoy平面上的不规则三角网；在不规则三角网上加上各点的高程信息后形成xyz空间不规则三角网(TIN)；

Delaunay三角剖分满足以下3个准则：

(1)最大化三角剖分T中所有三角形的最小角，以避免出现过于扁平的三角形；

(2)三角剖分T的任一三角形的外接圆为空圆，即该外接圆内不包含任何其他点，达到局部优化处理；

(3)四边形对应的剖分对角线两边的角度之和不大于180°。

步骤4：利用下式计算不规则三角网中每个空间三角形的形变量并保存；

$s=C\×\underset{i=1}{\overset{3}{\Σ}}|\sin A_i^{\′}-\sin A_i|---\left(2\right)$

其中，s是空间三角形的形变量，C是常数(C是为了增大s的取值范围，取大于1的任意值)，A_i表示xoy平面三角形的内角，A′_i表示添加高程信息形成空间三角形后与A_i对应的内角，i＝1,2,3。

步骤5：计算三角形的形变量阈值，具体步骤如下：

501：读取所有空间三角形的形变量s组成的形变量集合S；设定形变量阈值T^k的初始值T⁰：

$T^0=\frac{S_{\min }+S_{max}}{2}---\left(3\right)$

其中，S_min和S_max分别表示形变量集合S中的最小值和最大值；

502：将形变量集合S中所有的三角形形变量s分成形变量大于T^k的集合A和形变量小于等于T^k的集合B两部分；

503：利用式4、式5计算A、B两部分的平均形变量g_A和g_B：

$g_A=\frac{1}{n}\underset{S\left(i\right)>T^k}{\Σ}S\left(i\right)---\left(4\right)$

$g_B=\frac{1}{m}\underset{S\left(i\right)\≤T^k}{\Σ}S\left(i\right)---\left(5\right)$

其中，S(i)是形变量集合S中的第i个三角形形变量，n和m分别是A和B中三角形形变量s的个数；

504：利用式6计算新的阈值T^k+1：

$T^{k+1}=\frac{g_A+g_B}{2}---\left(6\right)$

505：判断条件|T^k+1-T^k|≤δ,k＝0,1,2,…是否成立，δ是允许误差，此处误差范围C*0.01-C*0.05，其中，C是式2中的常数；条件成立，则形变量阈值T^k更新为步骤504得到的T^k+1，然后执行步骤6；否则，阈值仍为T^k，此时令k＝k+1，并返回步骤503。

步骤6：将步骤4得到的形变量大于形变量阈值T^k的空间三角形标记为边缘三角形；如果形变量大于形变量阈值,进入步骤8，否则进入步骤7；

步骤7：计算狭长三角形：

利用式7计算步骤3得到的空间三角网中所有三角形的狭长度；比较各三角形的狭长度L_i和狭长度阈值Thresd，如果L_i大于Thresd，Thresd不小于100；则该三角形为狭长三角形，将该三角形标记为边缘三角形；

L_i＝|A_iG_i|²+|B_iG_i|²+|C_iG_i|² (7)

其中，L_i表示三角网中第i个三角形的狭长度；A_i，B_i，C_i分别为第i个三角形的三个顶点，G_i为第i个三角形的重心。

步骤8：确定边缘点，得到LiDAR点云数据的边缘图像，具体步骤如下：

801：粗边缘点的确定：

对于步骤6得到的每个边缘三角形利用式8计算高程阈值H(G)，并分别对每个边缘三角形中的顶点进行粗边缘点的选取：

$H\left(G\right)=\frac{H\left(A\right)+H\left(B\right)+H\left(C\right)}{3}---\left(8\right)$

其中，H(A)、H(B)、H(C)分别是A、B、C点的高程信息，单位：m；H(G)为高程阈值，单位：m；

若边缘三角形的三个顶点A,B,C中有任意两点A,B的高程均大于H(G)且|H(A)-H(B)|≤σ，σ取0.5-2.5，则将A,B两点均取为粗边缘点，否则，将A,B,C三点中高程最大的点取为粗边缘点。

对步骤7得到的每个边缘三角形的三个顶点都确定为粗边缘点。

:802：删除得到的所有粗边缘点中的孤立点，孤立点利用步骤2的公式9进行计算；将剩余的粗边缘点作为最终的边缘点输出，即得到LiDAR点云数据的边缘图像。

n＝∑(|P_i-P_j|＜σ) (9)

其中，P_i，P_j表示粗边缘点中点i和点j的坐标，n表示粗边缘点中与坐标P_i的距离小于σ的点的个数，n取0-10，σ取10-20m。

为了说明本发明的效果，发明人进行了如下仿真试验。

仿真试验1：

对本发明中无数据空白的LiDAR点云数据进行边缘检测仿真效果。仿真条件：MATLAB7.0软件。遵循上述本发明的技术方案，参见图2，对标准数据Sample23进行边缘检测仿真。本试验中，步骤2中n的取值为6，σ的取值为15m；步骤4和步骤5中的常数C的取值为1000；步骤5中δ的取值为10；步骤6中狭长三角形的阈值为100；步骤7中σ的取值为1.5m。

从图2(d)可以看出，本发明的方法保留了LiDAR点云的全部信息，有利于后续的特征提取、配准融合等处理。而将LiDAR点云数据进行深度图像化后利用对灰度图像进行边缘检测的方法会丢失一些数据，灰度化过程中会减少一些边缘信息，而且没有保持LiDAR点云数据的有效信息，不利于后续的处理。

仿真试验2：

对本发明中含有河流等数据空白区域的LiDAR数据进行边缘检测的仿真。仿真条件：MATLAB7.0软件。遵循上述本发明的技术方案，参照图3，对标准数据Sample61进行边缘检测仿真。本试验中，步骤2中n的取值为6，σ的取值为15m；步骤4和步骤5中常数C的取值为1000；步骤5中δ的取值为10；步骤6中狭长三角形的阈值为100；步骤7中σ的取值为1.5m。

从图3(d)可以看出，本发明的方法能够很好的检测到数据空白区域的边缘点，利用狭长三角形作为测度检测到河流边缘点。而将LiDAR点云数据进行深度图像化后进行插值会将空白区域消除，再利用灰度图像的边缘检测方法进行检测时会出现检测不到河流边缘或者检测到一些错误的边缘。

Claims (6)

1.一种基于不规则三角网的LiDAR点云数据边缘检测方法，其特征在于，具体包括如下步骤：

步骤1：读取目标区域原始LiDAR点云数据；

步骤2：剔除粗差噪声点；

步骤3：将步骤2得到的剔除粗差噪声后的LiDAR点云数据进行Delaunay三角剖分，生成xoy平面上的不规则三角网；在不规则三角网上加上各点的高程信息后形成xyz空间不规则三角网TIN；

步骤4：计算不规则三角网中每个空间三角形的形变量并保存；

步骤5：计算三角形的形变量阈值；

步骤6：将步骤4得到的形变量大于形变量阈值的空间三角形标记为边缘三角形；如果形变量大于形变量阈值,进入步骤8，否则进入步骤7；

步骤7：计算狭长三角形，比较各三角形的狭长度L_i和狭长度阈值Thresd，如果L_i大于Thresd，则该三角形为狭长三角形，将该三角形标记为边缘三角形，其中，Thresd不小于100；

步骤8：确定边缘点，得到LiDAR点云数据的边缘图像，具体步骤如下：

801：粗边缘点的确定：

对于步骤6得到的每个边缘三角形利用式1计算高程阈值H(G)，并分别对每个边缘三角形中的顶点进行粗边缘点的选取：

其中，H(A)、H(B)、H(C)分别是A、B、C点的高程信息，单位：m；H(G)为高程阈值，单位：m；

若边缘三角形的三个顶点A,B,C中有任意两点A,B的高程均大于H(G)且|H(A)-H(B)|≤σ，σ取0.5-2.5m，则将A,B两点均取为粗边缘点，否则，将A,B,C三点中高程最大的点取为粗边缘点；

对步骤7得到的每个边缘三角形的三个顶点都确定为粗边缘点；

802：删除所有粗边缘点中的孤立点，将剩余的粗边缘点作为最终的边缘点输出， 即得到LiDAR点云数据的边缘图像。

2.如权利要求1所述的基于不规则三角网的LiDAR点云数据边缘检测方法，其特征在于，所述步骤2的粗差噪声点利用式2计算得到：

n＝∑(|P_i-P_j|＜σ) (2)

其中，P_i表示原始LiDAR点云数据中点i的XY坐标，n表示原始LiDAR点云数据中与坐标P_i的距离小于σ的点的个数，n取0-10，σ取10-20m。

3.如权利要求1所述的基于不规则三角网的LiDAR点云数据边缘检测方法，其特征在于，所述步骤4的计算不规则三角网中每个空间三角形的形变量是指：对每个三角形利用式3进行计算：

其中，s是空间三角形的形变量，C是常数，取大于1的任意值，A_i表示xoy平面三角形的内角，A′_i表示添加高程信息形成空间三角形后与A_i对应的内角，i＝1,2,3。

4.如权利要求1所述的基于不规则三角网的LiDAR点云数据边缘检测方法，其特征在于，所述步骤5计算三角形的形变量阈值的具体步骤如下：

501：读取所有空间三角形的形变量s组成的形变量集合S；设定形变量阈值T^k的初始值T⁰：

其中，S_min和S_max分别表示形变量集合S中的最小值和最大值；

502：将形变量集合S中所有的三角形形变量s分成形变量大于T^k的集合A和形变量小于等于T^k的集合B两部分；

503：利用式5、式6计算A、B两部分的平均形变量g_A和g_B：

其中，S(i)是形变量集合S中的第i个三角形形变量，n和m分别是A和B中三角形形变量s的个数；

504：利用式7计算新的阈值T^k+1：

505：判断条件|T^k+1-T^k|≤δ,k＝0,1,2,…是否成立，δ是允许误差，此处误差范围C*0.01-C*0.05，其中，C是式2中的常数；条件成立，则形变量阈值T^k更新为步骤504得到的T^k+1，然后执行步骤6；否则，阈值仍为T^k，此时令k＝k+1，并返回步骤503。

5.如权利要求1所述的基于不规则三角网的LiDAR点云数据边缘检测方法，其特征在于，所述步骤7的计算狭长三角形是指：

利用式8计算步骤3得到的空间三角网中所有三角形的狭长度；

L_i＝|A_iG_i|²+|B_iG_i|²+|C_iG_i|² (8)

其中，L_i表示三角网中第i个三角形的狭长度；A_i，B_i，C_i分别为第i个三角形的三个顶点，G_i为第i个三角形的重心。

6.如权利要求1所述的基于不规则三角网的LiDAR点云数据边缘检测方法，其特征在于，所述步骤8的802中所述孤立点利用步骤2的公式9进行计算：

n＝∑(|P_i-P_j|＜σ) (9)

其中，P_i，P_j表示粗边缘点中点i和点j的坐标，n表示粗边缘点中与坐标P_i的距离小于σ的点的个数，n取0-10，σ取10-20m。