基于特征差异的机载LiDAR侵蚀沟点云分割方法

许 颖，岳东杰，孟德友

（1.河南财经政法大学资源与环境学院，郑州 450046；2.河海大学地球科学与工程学院，南京 211100）

1 研究背景

侵蚀沟作为一种常见的土壤侵蚀类型，侵蚀量大，速度快，在土壤侵蚀中占据着重要的位置［1］。随着对地观测技术和点云数据处理能力的发展，机载激光雷达（Light Detection and Ranging，LiDAR）的出现，为获取高时空分辨率的地球空间信息提供了全新的技术手段［2－3］。但是，点云数据量、特征差异和自然地物的复杂性严重制约着机载LiDAR技术数据后处理的发展。不像建筑物、树木、道路、桥梁、电力线［4］等有着明显地表特征的地物［5］，侵蚀沟嵌入在地表中，很难形成明显的表面特征，并且它们的变化范围、形成过程和外观不规则性更是增加了检测的难度［6］。

在LiDAR技术出现的初期，Betts等［7］借助由Li-DAR数据构建的数字高程模型（DEM）提取局部溪谷的形态信息。Perroy等［8］应用机载和地面激光雷达相结合的方法进行沟蚀评估，由于地面激光雷达扫描侵蚀沟底部和侧部有一定的限制，因此机载激光雷达技术更适合大范围侵蚀沟监测。Evans等［9］提出从高精度DEM中提取侵蚀沟，并在侵蚀沟边缘内插出其深度。Mason等［10］针对点云利用边界检测算子得到高梯度值，此方法在复杂地形中有一定的局限性。Rutzinger等［11］基于最大曲率分割算法提出对象驱动方法进行兴趣点检测，算法适于深且宽的轮廓特征，对纹理特征复杂的沟壑地区效果不佳。

以上基于DEM对提取类似侵蚀沟的微地貌特征还有一定的局限性，基于点云本身的操作会更符合地物地形特征。基于此，本文提出一种基于机载LiDAR点云的不同尺度下的表面特征差异侵蚀沟点云分割方法。考虑到侵蚀沟点云信息位于表面特征复杂山区，改进基于点云的边缘分割算法，基于法向量和曲率的变化关系来进行侵蚀沟点云分割。

2 特征差异分割方法

特征差异分割方法的思想是基于侵蚀沟处于地表发生突变的地方，选取不同半径，得到不同表面特征参数，设定一定阈值分割得到侵蚀沟点云。曲面的几何信息是指能够反映曲面局部特性的数据量［12］，其主要包括曲面的法向量、高斯曲率、平均曲率和曲面变化度等［13］。本文将法向量和曲率作为点云的表面特征参数。

2.1 点云特征分析

侵蚀沟的地貌形态一般表现为纵剖面与所在坡面基本一致，横剖面一般为U或V字形，因此通过切沟两边的沟沿线及U字形切沟沟底的两边或V字形切沟的沟底线可很好地反映切沟的三维形态［14］。由于沟壁斗力，沟谷深壑，因而侵蚀沟在水力、重力的双动力混合作用下，往往具有随机性特征。总结来看，侵蚀沟的点云特征主要表现以下几个方面：①从侵蚀沟形态上来看，其沟沿线点云表现出不规则、不对称的特征，沟壁点云密度分布不均匀，增加了侵蚀沟提取的难度，沟沿线和沟底点云高程差别明显；②从与其他地物差别来看，比如与建筑物、树木、道路、桥梁等地物表现出的点云特征进行比较，侵蚀沟表现出的点云嵌入在地形特征中，很难形成明显的表面特征，并且一般侵蚀沟发生在山区，地形起伏比较大，地物特征不规则；③侵蚀沟的动态变化更是增大了检测的难度。

2.2 法向量与曲率估计

法向量反映了曲面在某点处的几何特征，原始的散乱点云本身没有法向量，点云数据的法向量主要是通过点的邻域来确定［15］。因此，点的法向量近似为邻域表面的切平面法向量估计问题，对于点云中的每个点pi＝( xi，yi，zi)T，获取由距离半径确定的k个相邻点，每个点pi的协方差矩阵可以表示为

式中：k代表点pi邻域中所有点的个数；p－表示邻域中所有点的质心；λj是协方差矩阵的第j个特征值；v→—j是第j个特征向量。

用数学方法计算出的法向量解决不了它的符号问题，因此得到法向量的方向是任意的。为了解决这个问题，我们需要进行法向量一致性调整，选取视点vp，使其满足方程：

2.3 特征差异原理

特征差异处理为无序大量三维点云数据的处理提供了一种高效率、多尺度方法，以点集P中的每个点pi为圆心，选取 2个不同的半径r1，r2（r1＞r2），由上述方法会分别得到2个不同的法向量n1( pi，r1)，n2( pi，r2)和曲率值，如图1所示。表面特征差异定义为

2.4 基于特征差异的侵蚀沟提取

基于特征差异的侵蚀沟分割方法的主要工作是在特征差异的估算上。本文邻域点的确定采用固定距离法，其中距离的选择可以根据点云分辨率（点间距的大小）来灵活确定，一般情况下，领域点数量设定在10～20之间比较好。

基于特征差异的侵蚀沟提取具体的算法如下：

图1 半径不同产生的特征差异Fig.1 Characteristic difference produced by different radii

（1）选取半径r1得到邻域点，采用主成分分析方法拟合得到点云法向量和曲率值，并对法向量进行一致性调整。

（2）选取半径r2得到邻域点，采用主成分分析方法拟合得到点云法向量和曲率值，并对法向量进行一致性调整（r1＞r2）。

（3）根据法向量和曲率的不同，计算每个点的特征差异值。

（4）根据侵蚀沟的点云特征，设定误差阈值（0.25），提取差异特征大的值。

（5）基于距离（3 m）聚类的方法，剔除小块特征差异大的点，进而得到侵蚀沟点云。

因此，基于不同尺度下的表面特征差异侵蚀沟分割算法流程如图2所示。

图2 点云分割流程Fig.2 Flow chart of point cloud segmentation

3 试验及结果分析

3.1 数 据

为了测试本文算法的可行性，在VC＋＋环境下实现了该算法。试验数据为甘肃省嘉峪关市某区域，选取区域共358 977个点中的一部分，经过噪声处理和滤波后地面点数为11 038，点云平均密度为14个点／m2，点云数据的三维显示如图3所示。

图3 点云三维显示图Fig.3 Three-dimensional display of point cloud

3.2 试验过程

本文根据固定距离法选择邻域点，点云密度14个点／m2，点云密度比较高，半径分别取1 m和2 m选择邻域点。试验过程如下所述。

（1）对地面点云进行归一化处理，提高计算效率。

（2）对每一点云取邻域半径分别为1 m和2 m，计算法向量和曲率值，进而得到特征差异值，设定误差阈值为0.02，得到分割后结果，如图4（a）所示。

（3）设定聚类距离为2 m，得到侵蚀沟点云，如图 4（b）所示。

图4 侵蚀沟点云提取过程Fig.4 Process of extracting point cloud of gully

对图3截取的原始点云进行极大曲率特征点提取［14］，提取结果如图5（a）所示，于是特征点与图 4分割得到的侵蚀沟点云合并得到完整沟壑点云，如图 5（b）所示。

基于提取的侵蚀沟点云，对其进行定量分析，分别得到侵蚀沟的长度约为32 m，宽度约为2.5 m，如图6所示。侵蚀沟定量数据的获得，为以后侵蚀沟定量化评价和水土流失侵蚀量计算提供了一种新的方法。

图5 完整沟壑点云Fig.5 Integrated point cloud of gully

图6 侵蚀沟参数计算Fig.6 Calculation of parameters of gully

3.3 结果分析

视觉上来看，本文提取的侵蚀沟点云与原始点云形态相似，符合沟的形态特征，如图5所示。另一方面，对原始三维点云和分割得到的沟壑点云叠加，如图7（a）所示，可以看出所分割点云完全落在原始三维点云的沟壑位置上，符合沟壑点云特征；对原始点云建立DEM，通过侵蚀沟三维点云与DEM叠加，如图7（b）所示，可以看出，侵蚀沟沟沿线与DEM的高程变化相吻合。以上试验结果表明，采用特征差异方法能从点云中准确分割出比较完整的侵蚀沟点云信息。

图7 点云分割结果验证Fig.7 Verification of point cloud segmentation results

4 结 论

本文提出的基于特征差异的机载LiDAR侵蚀沟点云分割方法，考虑到侵蚀沟点云特征以及所处地形特性的特殊性，从仅仅考虑点云法向量和曲率单因素微分信息扩展到邻域的多尺度分析，即从它们变化邻域的尺寸可以得到更多的曲面信息。从结果分析来看，多尺度的算法可以从散乱点云中分割精确的侵蚀沟点云，很好地保持了侵蚀沟形态特征。所提取的侵蚀沟点云能合理表达真实形态，为进一步进行侵蚀量评价提供基础。

［1］ WU Y，CHENG H.Monitoring of Gully Erosion on the Loess Plateau of China Using a Global Positioning System［J］.Catena，2005，63（2／3）：154－166.

［2］ 郑辑涛，张 涛.基于可变半径圆环和 B样条拟合的机载 LiDAR点云滤波［J］.测绘学报，2015，44（12）：1359－1366.

［3］ 张 煜，窦延娟，张晓东.机载激光雷达数据采集及数据处理［J］.长江科学院院报，2010，27（1）：13－16.

［4］ 韩文军.机载LiDAR数据的电力塔建模方法研究［J］.长江科学院院报，2012，29（10）：122－126.

［5］ 张春亢，赵学胜，王洪斌.采用Morse理论的小尺度地形特征提取方法［J］.测绘科学技术学报，2015，（3）：266－270.

［6］ FILIN S，AVNI Y，BARUCH A，et al.Characterization of Land Degradation along the Receding Dead Sea Coastal Zone Using Airborne Laser Scanning［J］.Geomorphology，2014，206：403－420.

［7］ BETTSH D，DEROSE R C.Digital Elevation Models as a Tool for Monitoring and Measuring Gully Erosion［J］.International Journal of Applied Earth Observation and Geoinformation，1999，1（2）：91－101.

［8］ PERROY R L，BOOKHAGEN B，ASNER G P，et al.Comparison of Gully Erosion Estimates Using Airborne and Ground-based LiDAR on Santa Cruz Island，California［J］.Geomorphology，2010，118（3）：288－300.

［9］ EVANSM，LINDSAY J.High Resolution Quantification of Gully Erosion in Upland Peatlands at the Landscape Scale［J］.Earth Surface Processes and Landforms，2010，35（8）：876－886.

［10］MASON D C，SCOTT T R，WANG H J.Extraction of Tidal Channel Networks from Airborne Scanning Laser Altimetry［J］.ISPRS Journal of Photogrammetry＆Remote Sensing，2006，61（2）：67－83.

［11］RUTZINGER M，HÖFLE B，PFEIFER N，et al.Objectbased Analysis of Airborne Laser Scanning Data for Natural Hazard Purposes Using Open Source Components［C／CD］∥ISPRS.International Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences.Salzburg，Austria，July 4－5，2006：36（part 4／C42）.

［12］于海洋，余鹏磊，谢秋平，等.机载LiDAR数据建筑物顶面点云分割方法研究［J］.测绘通报，2014，（6）：20－23.

［13］刘 宇.基于微分信息的散乱点云拼合和分割［D］.武汉：华中科技大学，2008.

［14］BARUCH A，FILIN S.Detection of Gullies in Roughly Textured Terrain Using Airborne Laser Scanning Data［J］.ISPRSJournal of Photogrammetry and Remote Sensing，2011，66（5）：564－578.

［15］TANG C K，MEDIONI G.Curvature-augmented Tensor Voting for Shape Inference from Noisy 3D Data［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2002，24（6）：858－864.