CN103810489A - 基于不规则三角网的LiDAR点云数据水上桥梁提取方法 - Google Patents

Abstract

本发明公开了一种基于不规则三角网的LiDAR点云数据水上桥梁提取方法：读取目标区域的原始LiDAR点云数据；剔除粗差噪声；对LiDAR点云数据进行Delaunay三角剖分生成不规则三角网；计算所有三角形的狭长度；计算边缘点高程阈值，以高程阈值为分界点分离桥梁边缘点和河流边缘点，得到河流边缘点和桥梁边缘点高程变化曲线；对河流边缘点进一步分离，拟合河岸边缘曲线；对桥梁边缘点进一步分离，拟合桥梁边缘曲线；求解桥梁角点，对得到的所有桥梁角点进行标记，提取出水上桥梁。该方法充分利用水体对激光点的吸收特性，对LiDAR点云数据进行三角剖分确定水体，通过分离河流边缘点和桥梁边缘点，得到桥梁角点，能够高效准确提取位于水体上的一个或多个桥梁。

Description

基于不规则三角网的LiDAR点云数据水上桥梁提取方法

技术领域

本发明属于遥感应用研究领域，具体涉及一种基于不规则三角网的LiDAR点云数据水上桥梁提取方法。 

背景技术

机载LiDAR系统能够直接获取地面三维数据，具有高精度、高密度、高效率和成本低等优点，利用它快速获得的海量数据可以得到我们所需要的各种图像产品，因此在现代测绘中发挥越来越重要的作用。桥梁作为重要的人工建筑和交通枢纽，利用LiDAR点云数据进行桥梁目标提取具有重要意义。LiDAR点云数据是指机载LiDAR系统通过发射和接收激光脉冲得到地表的高精度3维点坐标。 

桥梁目标的提取基本上是基于影像数据，基于LiDAR点云数据的桥梁提取算法较少。Sithole等人在文献“Sithole G,Vosselman G.‘Bridge detection in airborne laser scanner data’[J].ISPRS journal of photogrammetry and remote sensing,2006,61(1):33-46.”中提出的基于剖面分析的桥梁点提取算法，该方法能够较好的提取陆地上的任意形状的桥梁点，但对于水上桥梁的提取由于要先进行水体的提取等工作从而降低了水上桥梁的提取效率。李云帆等人在文献“李云帆,马洪超,邬建伟,等.‘顾及地形特征的机载LiDAR数据桥梁提取算法研究’[J].武汉大学学报:信息科学版,2011,36(5):552-555.”中提出的算法能够较好的提取水上和陆上桥梁，但是由于频繁的进行三角形定位和插入、删除操作影响了识别的效率；另外，算法中使用过多的参数设置也会影响桥梁点提取的效果。 

综上，基于LiDAR点云数据的桥梁提取是一个较新的研究领域，且上述现有的桥梁提取技术存在识别和提取效率低，算法中使用过多参数影响提取效果等缺点。 

研究者发现，由于水体对激光点的吸收作用导致LiDAR点云数据会出现数据空白，在进行Delaunay三角剖分时会出现狭长的三角形，利用这些狭长三角形进行桥梁检测和提取桥梁的方法不失为一种创新。 

发明内容

针对现有的LiDAR点云数据水上桥梁提取方法的效率不高和参数设置较多等缺点，提出了基于不规则三角网LiDAR点云数据水上桥梁提取方法，该方法充分利用水体对激光点的吸收特性，对LiDAR点云数据进行三角剖分，在三角网空白区域检测狭长三 角形来确定水体，然后通过分离河流边缘点和桥梁边缘点，对边缘点进行边缘曲线拟合，得到桥梁角点，最终能够高效准确提取位于水体上的一个或多个桥梁，具有很好的实际应用。本发明的可以应用于LiDAR点云数据中有水桥梁的识别和提取，从而将获取到的桥梁信息用于交通方面以及后续的三维重建。 

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

一种基于不规则三角网的LiDAR点云数据水上桥梁提取方法，具体包括如下步骤： 

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

步骤2，剔除粗差噪声； 

步骤3，对步骤2得到的剔除粗差后的LiDAR点云数据进行Delaunay三角剖分：将离散分布的点云数据生成xoy平面上的不规则三角网； 

步骤4，计算步骤3得到的三角形网中所有三角形的狭长度，并提取边缘点。 

步骤5，计算边缘点高程阈值，以高程阈值为分界点分离桥梁边缘点和河流边缘点，得到河流边缘点和桥梁边缘点高程变化曲线； 

步骤6，对步骤5中得到的河流边缘点进一步分离，拟合河岸边缘曲线； 

步骤7，对桥梁边缘点进一步分离，拟合桥梁边缘曲线； 

步骤8，求解桥梁角点，具体步骤如下： 

对步骤6、步骤7得到的两测河岸边缘的拟合曲线方程以及每个桥梁的每一侧边缘点的拟合曲线方程联立求解，得到所有桥梁的角点； 

联立方程如下： 

$\left\{\begin{array}{cc}{\mathrm{RL}}_i(x,y)=0& i=\mathrm{1,2}\\ {\mathrm{BL}}_j(x,y)=0& j=1,.2\·\·m\end{array}\right.$

其中，RL_i(x,y)＝0,i＝1,2，指步骤6中拟合的两侧河岸的边缘点的拟合曲线，其中，i表示河流的两岸；BL_j(x,y)＝0,j＝1,2··m分别表示步骤7中拟合的不同桥梁的各边缘点的拟合曲线，其中，j表示不同的桥梁边缘，m表示提取的桥梁边缘的个数最大值。 

步骤9，对步骤8得到的所有桥梁角点进行标记，提取出水上桥梁。 

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

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

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

进一步的，所述步骤4的具体步骤如下： 

401：计算步骤3得到的不规则三角网中所有三角形的狭长度L_i： 

L_i＝|A_iG_i|² ₊|B_iG_i|² ₊|C_iG_i|²

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

402：比较各三角形的狭长度L_i和狭长度阈值Thresd，如果狭长度L_i大于狭长度阈值Thresd，则该三角形为狭长三角形，将该三角形的各顶点作为边缘点提取出来。 

进一步的，所述步骤5的具体步骤如下： 

501：读取步骤4提取的所有边缘点，将所有边缘点的按高程排序，得到河流边缘点和桥梁边缘点高程变化曲线；将该曲线上第一个高程值突变点作为河流边缘点的高程变化边界点，将第二个高程值突变点作为桥梁边缘点的高程变化边界点； 

502：利用下式计算高程阈值Z，并将步骤4得到的边缘点中，所有高程大于高程阈值的边缘点标记为桥梁边缘点，否则将其标记为河流边缘点； 

$Z=\frac{H\left(A\right)+H\left(B\right)}{2}$

其中，Z表示高程阈值，单位：m；H(A)表示A点处的高程，单位：m；H(B)表示B点处的高程，单位：m；A点表示河流边缘点的高程变化的边界点，B点表示桥梁边缘点的高程变化的边界点。 

进一步的，所述步骤6的具体步骤如下： 

601：对步骤5中标记为河流边缘点进行Delaunay三角剖分，构建平面不规则三角网，然后利用步骤401中的公式计算得到的三角网中所有三角形的狭长度，将狭长度大于狭长度阈值thresd的三角形标记为狭长三角形，并计算所有狭长三角形的重心坐标。狭长三角形重心坐标的集合如下： 

G＝{(x,y)|L＞thresd} 

其中，L是指步骤401得到的三角形的狭长度； 

对所有狭长三角形的重心坐标（x，y）进行曲线拟合，得到曲线RG(x,y)＝0，此 处拟合使用的是matlab提供的最小二乘法曲线拟合方法；以曲线RG(x,y)＝0为对比，将满足RG(x_i,y_i)＞0的河流边缘点(x_i,y_i)标识为河流a岸的边缘点，将满足RG(x_i,y_i)＜0的河流边缘点(x_i,y_i)标识为河流b岸的边缘点，实现河流两岸的河流边缘点的分离。 

602：分别对河流a岸、河流b岸的所有河流边缘点进行曲线拟合，得到两条拟合曲线方程RL_i(x,y)=0，i=1,2，i表示河流的不同两岸。 

进一步的，所述步骤7的具体步骤如下： 

701：对步骤5中标记的桥梁边缘点进行判别，分离出不同的桥梁的边缘点以及同一桥梁两侧边的边缘点；具体方法如下： 

(1)取U(k)＝{}，k＝1，点到直线距离阈值ε； 

(2)将桥梁边缘点集合Bri中的点从小到大排序，先按x坐标从小到大排序，如果其中相邻两点距离超过阈值θ，θ取1～5，则按y坐标从小到大重新排序； 

(3)取桥梁边缘点集合Bri中按（2）得到的排序的最小三个点拟合得到直线L(x,y)＝0，i＝1，将同一桥梁同一侧的边缘点与其他边缘点进行区分并归类，以便后续进行桥梁边缘拟合，具体如下： 

(3a)如果桥梁边缘点集合Bri中，有点(x_i,y_i)到直线L(x,y)＝0的距离d≤ε，则U(k)＝U(k)∪{(x_i,y_i)}，i＝i+1； 

(3b)如果桥梁边缘点集合Bri中，有点(x_i,y_i)到直线L(x,y)＝0的距离d＞ε，则i＝1，B＝B-U(k)，k＝k+1，转（4）； 

(4)如果B＝{}，则停止，否则转（2）； 

其中，U(k)是每次迭代完成后得到的各个桥梁边缘点的集合；k与i是迭代过程中递增的变量，k表示所有桥梁边缘点的迭代，i表示同一桥梁同一侧的桥梁边缘点的迭代；由于桥梁边缘是直线或者近似于直线，所以ε取值很小，取1～5； 

方法结束时，k的值作为迭代最终得到的桥梁边缘点总数，桥梁的个数为k/2，U(k)中存放的是上述方法最终得到的各个桥梁边缘的LiDAR点坐标。 

702：分别对步骤701得到的每个桥梁的每一侧的桥梁边缘点进行曲线拟合，得到 多个拟合曲线方程BL_j(x,y)=0，j=1,2…m，j表示不同的桥梁边缘，m是桥梁边缘的个数。 

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

第一，本发明在检测河流及桥梁边缘点时，应用Delaunay三角剖分建立不规则三角网，由于水体对LiDAR点的吸收特性会产生LiDAR点数据空白区域，空白区域生成的TIN三角形具有狭长性，可通过检测三角形的狭长度来确定河流区域从而检测出河流和桥梁边缘点，而此处检测出的狭长三角形在后续方法中还可利用其重心拟合曲线来分离河流两岸的河流边缘点，克服了现有水上桥梁提取方法中检测水体步骤繁复、效率不高等缺点，使得本发明中方法提取桥梁信息具有高效准确的特性。 

第二，在检测出河流和桥梁边缘点后，应用了高程差异以及曲线拟合来分离这些边缘点，河流和桥梁边缘点通过高程的差异进行分离，高程阈值通过统计边缘点高程通过突变点高程来确定，河岸两侧的河流边缘点通过拟合水体狭长三角形的重心来分离，而桥梁边缘点的分离通过初始点拟合的粗桥梁边缘直线，再根据桥梁边缘线周围的边缘点距离该直线的距离很小这一条件，不断拟合修正桥梁边缘线来确定。克服了现有技术中设置过多参数影响提取效果的问题，可在设置简单合理的阈值下，准确提取出水体上的一个或多个桥梁，方法具有较强的实用性。 

附图说明

图1为本发明的桥梁提取方法流程图。 

图2为本发明的仿真试验1中单桥梁提取效果图。其中，图2(a)和图2(c)是含有一个桥梁的原始LiDAR点云数据；图2(b)和图2(d)是检测识别出的桥梁并将其标记出的效果图。 

图3为本发明的仿真试验2中含有两个桥梁提取效果图。其中，图3(a)是含有两个桥梁的原始LiDAR点云数据；图3(b)是河流和桥梁边缘点检测效果图；图3(c)是检测识别出的桥梁并将其标记出的效果图。 

图4为本发明的仿真试验2中含有多个（三个或三个以上）桥梁提取效果图。其中，图4(a)是含有多个桥梁原始LiDAR点云数据；图4(b)是河流和桥梁边缘点检测效果图；图4(c)是检测识别出的桥梁并将其标记出的效果图。 

图5为河流边缘点和桥梁边缘点高程变化曲线。 

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

具体实施方式

参照图1，本发明的基于不规则三角网的LiDAR点云数据水上桥梁提取方法，具体包括如下步骤： 

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

步骤2，剔除粗差噪声； 

原始点云数据会产生一些粗差噪声，为了得到准确的点云数据，首先要对粗差噪声点进行剔除，粗差噪声点是指限定范围内没有其他点存在的孤立噪声点，其定义如下： 

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

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

步骤3，建立不规则三角网 

对步骤2得到的剔除粗差后的LiDAR点云数据进行Delaunay三角剖分：将离散分布的点云数据生成xoy平面上的不规则三角网。 

Delaunay三角剖分是指满足以下准则的三角剖分： 

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

（2）三角剖分T的任何一个三角形的外接圆满足空圆性质：该圆内不包含任何其他点，以达到局部优化处理。 

（3）四边形对应的剖分对角线两边的角度之和小于或者等于180°，则该剖分满足Delaunay特性。 

步骤4，计算步骤3得到的三角形网中所有三角形的狭长度，并提取边缘点。 

401：计算不规则三角网中所有三角形的狭长度L_i： 

L_i＝|A_iG_i|²+|B_iG_i|²+|C_iG_i|²

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

402：比较各三角形的狭长度L_i和狭长度阈值Thresd，如果狭长度L_i大于狭长度阈值Thresd，则该三角形为狭长三角形，将该三角形的各顶点作为边缘点提取出来； 

步骤5，计算边缘点高程阈值，以高程阈值为分界点分离桥梁边缘点和河流边缘点，得到河流边缘点和桥梁边缘点高程变化曲线。具体步骤如下： 

501：读取步骤4提取的所有边缘点，将所有边缘点的按高程排序，得到河流边缘点和桥梁边缘点高程变化曲线（参见图5）；将该曲线上第一个高程值突变点作为河流边缘点的高程变化边界点，将第二个高程值突变点作为桥梁边缘点的高程变化边界点； 

502：利用下式计算高程阈值Z，并将步骤4得到的边缘点中，所有高程大于高程阈值的边缘点标记为桥梁边缘点，否则将其标记为河流边缘点； 

$Z=\frac{H\left(A\right)+H\left(B\right)}{2}$

其中，Z表示高程阈值（单位：m）；H(A)表示A点处的高程（单位：m）；H(B)表示B点处的高程（单位：m）；A点表示河流边缘点的高程变化的边界点，B点表示桥梁边缘点的高程变化的边界点； 

H(A)和H(B)的计算过程如下： 

从河流边缘点和桥梁边缘点高程变化曲线中可看出，由于在A点和B点均有高程突变，A点和B点的梯度均会很大，利用下面两个式子得到A点和B点处的高程H(A)、H(B)，具体如下： 

$\▿H\left(k\right)=H(k+1)-H\left(k\right)$

$m=\arg \underset{k}{\max }(\▿H\left(k\right))$

其中，第一个式子用于迭代计算按高程排序的所有边缘点中，每两个相邻边缘点的高程之间的梯度值；第二个式子是对得到的梯度值由大到小进行排序，排序后取梯度值最大的和次大的，并将它们对应的序号k依次赋给m，m中的两个序号k对应的边缘点分别为点A和点B，点A和点B的高程分别为H（A）和H（B）。 

步骤6，对步骤5中得到的河流边缘点进一步分离，拟合河岸边缘曲线。 

601：对步骤5中标记为河流边缘点进行Delaunay三角剖分，构建平面不规则三角网，然后利用步骤401中的公式计算得到的三角网中所有三角形的狭长度，将狭长度大于狭长度阈值thresd的三角形标记为狭长三角形，并计算所有狭长三角形的重心坐标。狭长三角形重心坐标的集合如下： 

G＝{(x,y)|L＞thresd} 

其中，L是指步骤401得到的三角形的狭长度； 

对所有狭长三角形的重心坐标（x，y）进行曲线拟合，得到曲线RG(x,y)＝0，此处拟合使用的是matlab提供的最小二乘法曲线拟合方法；以曲线RG(x,y)＝0为对比，将满足RG(x_i,y_i)＞0的河流边缘点(x_i,y_i)标识为河流a岸的边缘点，将满足RG(x_i,y_i)＜0的河流边缘点(x_i,y_i)标识为河流b岸的边缘点，实现河流两岸的河流边缘点的分离。 

602：分别对河流a岸、河流b岸的所有河流边缘点进行曲线拟合，得到两条拟合曲线方程RL_i(x,y)=0，i=1,2（i表示河流的不同两岸），此处的拟合方法依然使用matlab提供的最小二乘法曲线拟合方法； 

步骤7，对桥梁边缘点进一步分离，拟合桥梁边缘曲线； 

701：对步骤5中标记的桥梁边缘点进行判别，分离出不同的桥梁的边缘点以及同一桥梁两侧边的边缘点；具体方法如下： 

(1)取U(k)＝{}，k＝1，点到直线距离阈值ε； 

(2)将桥梁边缘点集合Bri中的点从小到大排序，排序时先按x坐标从小到大排序，排序后，如果任一相邻两点距离

超过阈值θ，θ取1～5，证明两点的Y坐标相差很大，证明这样排序可能不合理，此时应按照y坐标从小到大重新排序； 

(3)取桥梁边缘点集合Bri中按（2）得到的排序的最小三个点拟合得到直线L(x,y)＝0，i＝1，将同一桥梁同一侧的边缘点与其他边缘点进行区分并归类，以便后续进行桥梁边缘拟合，具体如下： 

(3a)如果桥梁边缘点集合Bri中，有点(x_i,y_i)到直线L(x,y)＝0的距离d≤ε，则U(k)＝U(k)∪{(x_i,y_i)}，i＝i+1； 

(3b)如果桥梁边缘点集合Bri中，有点(x_i,y_i)到直线L(x,y)＝0的距离d＞ε，则i＝1，B＝B-U(k)，k＝k+1，转（4）； 

(4)如果B＝{}，则停止，否则转（2）； 

其中，U(k)是每次迭代完成后得到的各个桥梁边缘点的集合；k与i是迭代过程中递增的变量，k表示所有桥梁边缘点的迭代，i表示同一桥梁同一侧的桥梁边缘点的迭 代；ε取1～5； 

结束时，k的值作为迭代最终得到的桥梁边缘点总数，桥梁的个数为k/2，U(k)中存放的是上述方法最终得到的各个桥梁边缘的LiDAR点坐标； 

702：分别对步骤701得到的每个桥梁的每一侧的桥梁边缘点进行曲线拟合，得到多个拟合曲线方程BL_j(x,y)=0j=1,2…m（j表示不同的桥梁边缘，m是桥梁边缘的个数），此处的拟合方法依然使用的是matlab提供的最小二乘法曲线拟合方法； 

步骤8，求解桥梁角点。 

对步骤6、步骤7得到的两测河岸边缘的拟合曲线方程以及每个桥梁的每一侧边缘点的拟合曲线方程联立求解，得到所有桥梁的角点； 

联立方程如下： 

$\left\{\begin{array}{cc}{\mathrm{RL}}_i(x,y)=0& i=\mathrm{1,2}\\ {\mathrm{BL}}_j(x,y)=0& j=1,.2\·\·m\end{array}\right.$

其中，RL_i(x,y)＝0,i＝1,2，指步骤6中拟合的两侧河岸的边缘点的拟合曲线，其中，i表示河流的两岸；BL_j(x,y)＝0,j＝1,2··m分别表示步骤7中拟合的不同桥梁的各边缘点的拟合曲线，其中，j表示不同的桥梁边缘，m表示提取的桥梁边缘的个数最大值。 

步骤9，对步骤8得到的所有桥梁角点进行标记，提取出水上桥梁。 

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

仿真试验1： 

应用本发明的方法，对含有一个桥梁的LiDAR点云数据进行水上桥梁提取。仿真条件：MATLAB7.0软件。参照图2，对MARS7中样本LiDAR点云数据截取进行仿真实验。 

首先读入目标区域的原始LiDAR点云数据；对点云数据进行粗差剔除，该步骤中阈值n=3，σ=15m；之后对点云数据进行三角剖分，建立不规则三角网；对所生成三角网中的三角形计算三角形狭长度，提取狭长三角形各顶点作为边缘点，该步骤中狭长度阈值Thresd=100；将边缘点中的河流边缘点和桥梁边缘点进行分离（高程阈值Z如步骤5中计算得出，为30）；将河流边缘点进一步分离，拟合河岸边缘曲线该步骤中计算河流边缘点生成的狭长三角形时，狭长三角形阈值thresd设为100；将桥梁边缘点进 一步分离，拟合桥梁曲线，该步骤中距离阈值ε=3，相邻两点距离超过阈值θ=3；联立拟合出的河流边缘曲线方程和桥梁边缘曲线方程，求解所有桥梁角点；最终对提桥梁角点进行标记，提取水上桥梁。 

从图2(b)和图2(d)可以看出，本发明的方法效率高用时短，提取的桥梁角点误差小，提取效果好，能够高效准确地从含有单个桥梁的LiDAR点云数据中提取出桥梁。 

仿真试验2： 

应用本发明的方法，对含有两个桥梁的LiDAR点云数据进行水上桥梁提取。仿真条件：MATLAB7.0软件。参照图3，对MARS7中样本LiDAR点云数据截取进行仿真实验。 

首先读入目标区域的原始LiDAR点云数据；对点云数据进行粗差剔除，该步骤中阈值n=0，σ=15m；之后对点云数据进行三角剖分，建立不规则三角网；对所生成三角网中的三角形计算三角形狭长度，提取狭长三角形各顶点作为边缘点，该步骤中狭长度阈值Thresd=100；将边缘点中的河流边缘点和桥梁边缘点进行分离（高程阈值Z如步骤5中计算得出，为14）；将河流边缘点进一步分离，拟合河岸边缘曲线（该步骤中计算河流边缘点生成的狭长三角形时，狭长三角形阈值thresd设为100）；将桥梁边缘点进一步分离，拟合桥梁曲线（该步骤中距离阈值ε=3，相邻两点距离超过阈值θ=3）；联立拟合出的河流边缘曲线方程和桥梁边缘曲线方程，求解所有桥梁角点；最终对提桥梁角点进行标记，提取水上桥梁。 

从图3(c)可以看出，本发明的方法效率高用时短，提取的桥梁角点误差小，提取效果好，能够高效准确的从只含有两个桥梁的LiDAR点云数据中提取出桥梁。 

仿真试验3： 

应用本发明的方法，对含有三个桥梁的LiDAR点云数据进行水上桥梁提取。仿真条件：MATLAB7.0软件。参照图4，对MARS7中样本LiDAR点云数据截取进行仿真实验。 

首先读入目标区域的原始LiDAR点云数据；对点云数据进行粗差剔除，该步骤中阈值n=0，σ=15m；之后对点云数据进行三角剖分，建立不规则三角网；对所生成三角网中的三角形计算三角形狭长度，提取狭长三角形各顶点作为边缘点，该步骤中狭长度阈值Thresd=100；将边缘点中的河流边缘点和桥梁边缘点进行分离（高程阈值Z如步骤5中计算得出，为16）；将河流边缘点进一步分离，拟合河岸边缘曲线，该步骤中 计算河流边缘点生成的狭长三角形时，狭长三角形阈值thresd设为100；将桥梁边缘点进一步分离，拟合桥梁曲线，该步骤中距离阈值ε=3，相邻两点距离超过阈值θ=3；联立拟合出的河流边缘曲线方程和桥梁边缘曲线方程，求解所有桥梁角点；最终对提桥梁角点进行标记，提取水上桥梁。 

从图4(c)可以看出，本发明的方法效率高用时短，提取的桥梁角点误差小，提取效果好，能够高效准确的从只含有三个桥梁的LiDAR点云数据中提取出桥梁。 

Claims (6)

1.一种基于不规则三角网的LiDAR点云数据水上桥梁提取方法，其特征在于，具体包括如下步骤： 

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

步骤2，剔除粗差噪声； 

步骤3，对步骤2得到的剔除粗差后的LiDAR点云数据进行Delaunay三角剖分：将离散分布的点云数据生成xoy平面上的不规则三角网； 

步骤4，计算步骤3得到的三角形网中所有三角形的狭长度，并提取边缘点。 

步骤5，计算边缘点高程阈值，以高程阈值为分界点分离桥梁边缘点和河流边缘点，得到河流边缘点和桥梁边缘点高程变化曲线； 

步骤6，对步骤5中得到的河流边缘点进一步分离，拟合河岸边缘曲线； 

步骤7，对桥梁边缘点进一步分离，拟合桥梁边缘曲线； 

步骤8，求解桥梁角点，具体步骤如下： 

对步骤6、步骤7得到的两测河岸边缘的拟合曲线方程以及每个桥梁的每一侧边缘点的拟合曲线方程联立求解，得到所有桥梁的角点； 

联立方程如下： 

其中，RL_i(x,y)＝0,i＝1,2，指步骤6中拟合的两侧河岸的边缘点的拟合曲线，其中，i表示河流的两岸；BL_j(x,y)＝0,j＝1,2··m分别表示步骤7中拟合的不同桥梁的各边缘点的拟合曲线，其中，j表示不同的桥梁边缘，m表示提取的桥梁边缘的个数最大值。 

步骤9，对步骤8得到的所有桥梁角点进行标记，提取出水上桥梁。 

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

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

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

3.如权利要求1所述的基于不规则三角网的LiDAR点云数据水上桥梁提取方法，其特征在于，所述步骤4的具体步骤如下： 

401：计算步骤3得到的不规则三角网中所有三角形的狭长度L_i： 

L_i＝|A_iG_i|²+|B_iG_i|²+|C_iG_i|²

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

402：比较各三角形的狭长度L_i和狭长度阈值Thresd，如果狭长度L_i大于狭长度阈值Thresd，则该三角形为狭长三角形，将该三角形的各顶点作为边缘点提取出来。 

4.如权利要求1所述的基于不规则三角网的LiDAR点云数据水上桥梁提取方法，其特征在于，所述步骤5的具体步骤如下： 

501：读取步骤4提取的所有边缘点，将所有边缘点的按高程排序，得到河流边缘点和桥梁边缘点高程变化曲线；将该曲线上第一个高程值突变点作为河流边缘点的高程变化边界点，将第二个高程值突变点作为桥梁边缘点的高程变化边界点； 

502：利用下式计算高程阈值Z，并将步骤4得到的边缘点中，所有高程大于高程阈值的边缘点标记为桥梁边缘点，否则将其标记为河流边缘点； 

其中，Z表示高程阈值，单位：m；H(A)表示A点处的高程，单位：m；H(B)表示B点处的高程，单位：m；A点表示河流边缘点的高程变化的边界点，B点表示桥梁边缘点的高程变化的边界点。 

5.如权利要求1所述的基于不规则三角网的LiDAR点云数据水上桥梁提取方法，其特征在于，所述步骤6的具体步骤如下： 

601：对步骤5中标记为河流边缘点进行Delaunay三角剖分，构建平面不规则三角网，然后利用步骤401中的公式计算得到的三角网中所有三角形的狭长度，将狭长度大于狭长度阈值thresd的三角形标记为狭长三角形，并计算所有狭长三角形的重心坐标。狭长三角形重心坐标的集合如下： 

G＝{(x,y)|L＞thresd} 

其中，L是指步骤401得到的三角形的狭长度； 

对所有狭长三角形的重心坐标（x，y）进行曲线拟合，得到曲线RG(x,y)＝0，此 处拟合使用的是matlab提供的最小二乘法曲线拟合方法；以曲线RG(x,y)＝0为对比，将满足RG(x_i,y_i)＞0的河流边缘点(x_i,y_i)标识为河流a岸的边缘点，将满足RG(x_i,y_i)＜0的河流边缘点(xi,yi)标识为河流b岸的边缘点，实现河流两岸的河流边缘点的分离。 

602：分别对河流a岸、河流b岸的所有河流边缘点进行曲线拟合，得到两条拟合曲线方程RL_i(x,y)=0，i=1,2，i表示河流的不同两岸。 

6.如权利要求1所述的基于不规则三角网的LiDAR点云数据水上桥梁提取方法，其特征在于，所述步骤7的具体步骤如下： 

701：对步骤5中标记的桥梁边缘点进行判别，分离出不同的桥梁的边缘点以及同一桥梁两侧边的边缘点；具体方法如下： 

(1)取U(k)＝{}，k＝1，点到直线距离阈值ε； 

(2)将桥梁边缘点集合Bri中的点从小到大排序，先按x坐标从小到大排序，如果其中相邻两点距离超过阈值θ，θ取1～5，则按y坐标从小到大重新排序； 

(3)取桥梁边缘点集合Bri中按（2）得到的排序的最小三个点拟合得到直线L(x,y)＝0，i＝1，将同一桥梁同一侧的边缘点与其他边缘点进行区分并归类，以便后续进行桥梁边缘拟合，具体如下： 

(3a)如果桥梁边缘点集合Bri中，有点(x_i,y_i)到直线L(x,y)＝0的距离d≤ε，则U(k)＝U(k)∪{(x_i,y_i)}，i＝i+1； 

(3b)如果桥梁边缘点集合Bri中，有点(x_i,y_i)到直线L(x,y)＝0的距离d＞ε，则i＝1，B＝B-U(k)，k＝k+1，转（4）； 

(4)如果B＝{}，则停止，否则转（2）； 

其中，U(k)是每次迭代完成后得到的各个桥梁边缘点的集合；k与i是迭代过程中递增的变量，k表示所有桥梁边缘点的迭代，i表示同一桥梁同一侧的桥梁边缘点的迭代；ε取1～5； 

结束时，k的值作为迭代最终得到的桥梁边缘点总数，桥梁的个数为k/2，U(k)中存放的是上述方法最终得到的各个桥梁边缘的LiDAR点坐标； 

702：分别对步骤701得到的每个桥梁的每一侧的桥梁边缘点进行曲线拟合，得到多个拟合曲线方程BL_j(x,y)=0，j=1,2…m，j表示不同的桥梁边缘，m是桥梁边缘的个数。 

Images