基于无人机LiDAR数据的陡崖反坡滤波技术及应用

万年锋 袁启伦 邸国辉 周胜洁

摘要：

为了解决传统的不规则三角网（Triangular irregular network，TIN）算法在陡崖反坡地形滤波时，由于陡崖反坡会产生多个投影交点，坡度出现负值，陡崖反坡区域岩壁点往往被错误过滤的问题，提出一种刚性坐标变换与TIN结合的机载激光雷达（Light detection and ranging，LiDAR）点云滤波算法。采用刚性坐标变换的方法重新构建陡崖面的直观形态，将陡崖反坡等特异地形常规化。以湖北省恩施市姚家平水利枢纽工程坝址区域的陡崖反坡数据进行滤波分析，构建了侧向数字高程模型（DEM）表达陡崖地形，实现了侧向DEM在坝体填筑量计算上的应用。结果表明：该算法能够在陡崖反坡区域大幅降低Ⅰ类误差，并将Ⅱ类误差控制在一定范围内；应用该技术可为大坝设计提供直观准确的DEM，为三维设计提供了新思路。

关键词：

陡崖反坡地形； 点云滤波； 机载LiDAR； 不规则三角网； 侧向DEM； 大坝设计

中图法分类号：P237

文献标志码：A

DOI：10.15974/j.cnki.slsdkb.2023.05.007

文章编号：1006-0081（2023）05-0045-05

0 引 言

机载激光雷达（Light detection and ranging，LiDAR）技术结合了全球定位系统、激光扫描仪与惯导系统，能够快速高效地获取高精度地面点云数据［1］，具有效率高、受影响小、精度高、风险小等优势［2-4］。由于原始点云是由不同地形特征（地面、建筑物等）返回的大量离散点组合而成，并无语义信息或拓扑关联，因此在应用之前，必须将地面和非地面点分离，这个过程称为点云滤波。Lindenberger［5］较早应用数学形态学对剖面式测高数据执行开运算，然后使用自回归过程调整运算结果，实现对原始数据滤波，该算法原理简单、实现效率高，但缺乏适应性，容易削平地形。Axelsson［6］提出了三角网加密滤波算法，采用不规则三角形拟合地形，不断将满足地形特征的点加入三角网，通过三角网加密的形式还原地形信息，该算法适应性较强，已被应用于滤波处理软件TerraSolid模块中，但易丢失地面点且低矮植被不能有效剔除。Kraus等［7］提出了最小二乘内插滤波算法，根据待判点与拟合曲面之间的高程差确定相应的权重，然后依据权重值对拟合曲面进行调整，反复以上操作，直到完成所有数据，该算法滤波精度较高，但需要多次迭代，运行速度慢。罗伊萍等［8］采用多尺度的方式有效改善了固定窗口大小，提升了数学形态学算法的适应性，但组织形式有损原始数据的精度。徐国杰等［9］采用设定点到三角面的最大距离和点与三角面顶点的连线形成的最大角的方式对不规则三角网进行改进，并对参数动态调整，多次迭代完成滤波，但运行效率较差。张皓等［10］使用虚拟网格组织数据，提出了一种改进坡度的滤波算法，可以有效避免降低數据精度，但阈值的自适应性较差。上述各种主流的LiDAR滤波与分类算法研究都有其各自特点及局限性。其中，三角网滤波算法在山区等坡度变化大的地区适应性最佳，总误差值相对较小，滤波时间也更短［11］，但是对于山区里的一些特殊地形，比如水利工程中常见的陡崖及反崖地形，运用三角网滤波算法时，会出现一些陡崖崖壁及反崖反向坡度的数据损失，无法真实表达其地理地貌［12］。本文首次提出一种刚性坐标变换与不规则三角网（triangular irregular network，TIN）相结合的机载LiDAR点云滤波算法，能在陡崖及反坡等特殊地形高效准确进行滤波的方法，并结合湖北省恩施市姚家平水利枢纽工程，对无人机载LiDAR在陡崖反坡地形测绘及大坝设计中的应用进行实际分析。

1 陡崖反坡滤波方法

1.1 陡崖反坡

目前LiDAR点云滤波软件均未考虑陡崖反坡存在的多个投影交点（如图1所示的 O，P，Q ）、坡度出现负值的情况，在点云滤波提取地面点时，陡崖区域岩壁点（ P和Q ）往往被忽略。本文提出了对陡崖反坡点云先进行变换处理再进行滤波的方法。

1.2 陡崖反坡点云变换处理

在扫描获取陡崖点云数据后，采用高精度定位定姿系统（Position and orientation system，POS）数据解算和激光点云数据解算法对获取的数据进行解析，得到工程坐标系（ XYH ）0的点云，并构建陡崖区域，剪裁后，得到陡崖区域点云文件LAS0。工程坐标系（ XYH ）0为项目工程指定的坐标系，该坐标系包 括X轴、Y轴和H轴 。

基于已有数据绘制出整个陡崖的上边缘线，并对上边缘线进行直线拟合，得到线段 AB，作法线垂直于AB。AB的 法线方向如图2所示。对点云LAS0进行的刚性坐标变换，包括第一次刚性变换和第二次刚性变换。

（1） 第一次刚性变换。对点云LAS0进行第一次刚性变换，将其纵坐标轴旋转 θ角至线段AB 的法线方向，与法线方向重合（图2），得到坐标系（ XYH ）1，其中坐标按式（1）换算：

x ′ y′ h′= cos θ－ sin θ0  sin θ cos θ0 001x y h（1）

式中： x y h为初始点云坐标；x′ y′ h′为第一次刚性变换后点云坐标；θ为 旋转角。

（2） 第二次刚性变换。将第一次刚性变换后的高程坐标作为新的 X轴坐标，第一次刚性变换后的X轴坐标转换为新的Y轴坐标，第一次刚性变换后的Y轴坐标转换为新的高程坐标，得到新坐标系下陡崖面的初始点云LAS1。以p为原始点云上的任一点，p′为第一次刚性变换后p点的对应点，p″为第二次刚性变换后p点的对应点。p点坐标为（x，y，h），p′坐标为（x′，y′，h′），p″坐标为（x″，y″，h″）。p绕原点旋转θ，得到点p′，相应的x，y坐标也旋转θ为x′，y′，高程h′等于h。再对p′进行第二次刚性变换：高程h′坐标值转换为第二次刚性变换后的X轴，即X″；X′坐标值转换为第二次刚性变换后的Y轴，即Y″；Y′坐标值转换为第二次刚性变换后的H轴，即H″，得到坐标系（XYH） 2的点云，即陡崖面的初始点云。坐标换算如式（2）所示：

x″ y″ h″=001 100 010x′ y′ h′（2）

第二次刚性变换后， x″的数值等于原高程h，而x″ 作为横坐标，投影方向与法线方向一致，消除了陡崖反坡形态，显著减小了陡崖面坡度，得到陡崖面的初始点云LAS1，可直观反映陡崖的形态。此外，由于消除了反坡形态，深度坐标面的坡度显著减小，实现了陡崖区域岩壁地表信息的高效提取。至此，LAS1的反坡被消除。

1.3 滤 波

利用点云处理软件TerraSolid，基于不规则三角网TIN加密滤波算法，先通过一些较低种子点生成一个稀疏TIN，然后考察每个点与TIN的距离和坡度，并逐层迭代加密，进行地面点和非地面点（植被）分类，由LAS1提取出地面点点云LAS2。

2 实例分析

2.1 研究区域概况

研究区位于湖北省恩施市姚家平水利枢纽，试验航摄区域面积约1 km2，测区平均海拔700 m，最大相对高差为800 m。测区沿河谷方向分布多处陡崖及反崖。

2.2 无人机LiDAR航飞

航飞采用华测BB4大黄蜂四旋翼无人机+AS-900HL多平台激光雷达系统，精度高、测程长，配备高精度光纤惯导系统，满足此次LiDAR扫描任务。

由于测区属于峡谷地区，高差大，故整个测区采用仿地飞行。仿地飞行需要的测区初始数字地表模型（Digital surface model，DSM）来源于前期0.2 m分辨率的倾斜摄影成果。为扫描陡崖和反坡区域的地形，在平行峡谷方向布设了4条航线。航线高度120 m，航向重叠率60%，旁向重叠均值35%，点云密度不少于8个/m2。测区初始DSM及仿地飞行航线如图3所示。

2.3 陡崖及反崖滤波

沿北岸峡谷方向的谷底有1条200 m长的陡崖反坡，其原始剖面點云见图4。

（1） 基于1.2节的刚性变换，得到新坐标系下陡崖面的初始点云LAS1，如图5所示。可明显观察到，变换处理消除了陡崖反坡形态，显著减小了陡崖面坡度。

（2） 利用点云处理软件TerraSolid，基于不规则三角网TIN加密滤波算法，进行地面点和非地面点（植被）分类，由LAS1提取出地面点点云LAS2。其中，反坡设置参数：坡度80°，建筑物长度10 m，反复角设置为12°，提取地面点。

至此，得到了经过空间转换后的滤波点云。若需要得到原坐标系下滤波点云，经过反向变换即可。此外，由于陡崖（特别是反崖）的特殊性，在构建DEM时仍会出现因为DEM数据水平面上某一点只能对应唯一的高程值而无法生成反坡处DEM的问题，导致陡崖反坡的真实地形无法表达。为此，提出构建侧向DEM的方法表达陡崖地形，更能满足实际工程设计的需要。

（3） 侧向DEM是利用Globalmapper软件，导入陡崖面地面点的点云LAS2生成的DEM，如图6所示。因为点云LAS2是经过变换旋转后的点云，数学关系如1.2节的点云变换原理所述，故称为侧向DEM。这种经过变换旋转后的DEM可以更好地表达反坡处地形，有利于后期设计。

2.4 结果分析

根据2003年国际摄影测量和遥感学会（ISPRS）发布的误差评价标准，对本次实验滤波结果进行定量评价［13］。此种误差评价方法根据激光点的误分率，大致可以分为3类：① Ⅰ类误差，即将地面点分类为非地面点的误差（拒真误差，omission error）；② Ⅱ类误差，即将非地面点分类为地面点的误差（纳伪误差，commission error）；③ Ⅲ类误差，即总误差。前两类反映了算法的适应性，总误差反映了算法的可行性。为了最大限度得到地形表面数据，应尽可能避免Ⅱ类误差［14］。

为便于对比分析，同时采用传统TIN算法及本文算法对陡崖反坡点云数据进行滤波实验和精度评价，结果见表1。与传统TIN算法相比，本文算法的Ⅰ类误差降低幅度最大，降低34.02％；Ⅱ类误差值大致不变；总误差降低16.65％。综上所述，本文算法在Ⅱ类误差无较大变化的同时，显著降低Ⅰ类误差，验证了该算法具有较好的普适性与有效性，能够获得更接近真实地形的高精度数字高程模型。

2.5 重力坝填筑量计算

结合恩施市姚家平水利枢纽工程的重力坝方案，采用侧向DEM计算坝体填筑量［15］。取LAS0的纵轴至坝体轴线方向，将坝体顺时针旋转90°，将坝体横断面导入侧向DEM，给定基点高程（左或右坝段长度），进行填筑量计算。

坝高97 m，内坡比1∶0.25，外坡比1∶0.15，底部宽度53.3 m，顶部宽度14.5 m，取基点高程425 m（相应左岸坝段长60 m），左岸坝段TIN地面点有8 661个，点密度为2.6个/m2。采用DEM棱柱法计算方量，得填筑量为210 489.68 m3，见图7。同理，右岸坝段也可通过这种方法计算，进而得到整个坝段的填筑量。

通常断面法土方量的计算采用由两相邻横断面的平均面积乘以两相邻横断面间距的方法（断面法土方计算平均公式），即由侧向DEM截取11个断面计算断面面积，计算所得填筑量为209 331.98 m3，与DEM棱柱法对比，相对误差0.55%。结果表明：DEM棱柱法计算方量的误差远小于SL 197-97《水电水利工程施工测量规范》的要求。

DEM棱柱法充分利用了陡崖面地面点的点云数据，计算方法严密，比断面法计算精度高，为三维设计提供了新方法。

3 结 论

（1） 本文在TIN算法基础上，提出了通过对点云进行空间变换，改变陡崖及反坡的直观形态的方法，改善了滤除效果和地形起伏区域的适应性，显著提高了滤波精度。与传统TIN算法相比，该算法能有效过滤出山区陡崖及反坡的真实地面点。

（2） 利用恩施市姚家平水利枢纽工程区域陡崖反坡机载LiDAR点云数据进行实例分析，其试验结果定量证明了刚性坐标变换与TIN结合的机载LiDAR点云滤波算法的精度及实用性。

（3） 在大坝设计方面，应用该方法生成侧向DEM进行重力坝填筑量计算，结果表明：DEM棱柱法计算方量的误差小于规范要求，可用于大坝三维设计。

参考文献：

［1］ 易志朝.机载LiDAR技术在水利工程大比例尺地形图测绘中的应用［J］.陕西水利，2022（5）：181-183.

［2］ 谈语，赵云毅，富博.复杂山区河道地形机载Lidar数据质量控制［J］.水利水电快报，2021，42（10）：22-26.

［3］ 刘世振，邓建华，冯国正，等.机载LiDAR在山区型河道地形测绘中的适用性研究［J］.人民长江，2021，52（1）：108-113.

［4］ 吴昊，张晓萌.机载LiDAR在长江中游河道测绘中的应用［J］.人民长江，2022，53（9）：109-114.

［5］ LINDENBERGER J.Laser-Profilmessungen zur topographischen Gelndeaufnahme［D］.Stuttgart：Stuttgart University，1993.

［6］ AXELSSON P.DEM generation from laser scanner data using adaptive TIN models［J］.International Archives of Photogrammetry & Remote Sensing，2000，23（B4）：110-117.

［7］ KRAUS K，PFEIFER N.Determination of terrain models in wooded areas with airborne laser scanner data［J］.ISPRS Journal of Photogrammetry & Remote Sensing，1998，53（4）：193-203.

［8］ 罗伊萍，姜挺，龚志辉，等.基于自適应和多尺度数学形态学的点云数据滤波方法 ［J］.测绘科学技术学报，2009，26（6）：426-429.

［9］ 徐国杰，胡文涛.一种改进的 LiDAR 点云 TIN 迭代滤波算法［J］.测绘地理信息，2010，35（1）：33-35.

［10］ 张皓，贾新梅，张永生，等.基于虚拟网格与改进坡度滤波算法的机载LIDAR数据滤波［J］.测绘科学技术学报，2009，26（3）：224-227，231.

［11］ 邹正，邹进贵，胡海洋.不同机载LiDAR点云滤波算法对比分析［J］.测绘地理信息，2021，46（5）：52-56.

［12］ 朱依民，田林亚，毕继鑫，等.基于PTD和改进曲面拟合的高山区水电工程机载激光雷达点云滤波方法［J］.水利水电科技进展，2021，41（1）：35-40.

［13］ 吉雨田，张春亢，尹耀.机载LiDAR点云自适应滤波算法［J］.测绘科学技术学报，2021，38（2）：142-147.

［14］ 马超，蒋好忱，秦先锋，等.LiDAR-DP点云自动分类及精度评定［J］.测绘与空间地理信息，2022，45（4）：240-242，246.

［15］ 孟永东，陈俊华.约束TIN在土石方工程三维辅助设计中的应用［J］.计算机应用，2010，30（增1）285-288.

（编辑：江 焘，高小雲）

Abstract：

In order to solve the problem that the steep cliff reverse slope will produce multiple projection intersections，the slope will have negative values，and the rock wall points in the steep cliff reverse slope area are often wrongly filtered when the traditional triangular irregular network （TIN） algorithm is used to filter the terrain of the steep cliff reverse slope，a light detection and ranging（LiDAR） point cloud filtering algorithm combined with rigid coordinate transformation and TIN was proposed.The intuitive shape of the steep cliff surface was reconstructed by rigid coordinate transformation，which could normalize the specific terrain such as the steep cliff reverse slope to facilitate further filtering.Lateral Digital Elevation Model（DEM） was constructed to express the steep cliff terrain，the application of lateral DEM on dam filling amount calculation was realized.The results showed that this algorithm could greatly reduce the class Ⅰ error in the steep cliff reverse slope area and control the class Ⅱ error within a certain range.This method could provide intuitive and accurate DEM for dam design，which provided a new method for three-dimensional design.

Key words：

steep cliff reverse slope； point cloud filtering； UAV LiDAR； triangular irregular network； lateral DEM； dam design