机载LiDAR点云中静态水体边界提取

周维娜，程晓光，严明，卢浩男

(1.常州市自然资源和规划服务中心，江苏 常州 213022；2.飞燕航空遥感技术有限公司，南京 210001)

0 引言

机载激光雷达(light detection and ranging，LiDAR)已经成为陆地地形测绘的主要手段之一。在机载LiDAR测绘内业生产中，水体断裂线是主要的断裂线类型之一[1-2]。可以简单认为，水体断裂线是水体在二维平面上的边界。除数字高程模型(digital elevation model，DEM)生产外，水体边界还可用于地图制图、河湖蓄水监测等领域。对于水体较多较小的测区，手工勾绘水体断裂线繁琐，工作量大。如何自动、准确地提取水体边界用作断裂线，对提高机载LiDAR内业生产效率很有意义。

用于陆地测绘的机载LiDAR一般采用波长在近红外波段的激光脉冲[3]。由于纯净水体在近红外波段反射率极低[4]，一般后向散射回波功率很弱，难以被激光接收器识别为信号，因而在点云中，水体一般呈现为空洞或点分布特别稀疏的情况。这使得基于机载LiDAR点云提取水体边界具备独特的优势。但是，由于风、水生植物、水质等因素的影响，水体内部经常存在一定数目的点，可能分布较密集、回波强度较高，甚至与岸边相连。

目前基于机载LiDAR点云提取水体边界的研究一般基于水面点稀少、回波强度弱等假设[5-7]，导致难以处理水体内点较多或回波较强的情况。有研究人员使用光学影像辅助水体边界提取[8-10]。但是在航摄光学影像上，不同水质、水深、水草和阴影的水体的特征不一致，且相邻图幅间可能存在色差，而卫星遥感影像的分辨率一般偏低。此外，提取建筑物屋顶边界等的方法[11]不适用于提取水体边界。

本文提出一种基于Delaunay三角网提取静态水体边界的方法。该方法不依赖影像辅助、点云预分类和回波强度假设，不需要已知水面高程，只需要水体存在面积够大的平面空洞，既适用于水面点较少的情况，也适用于水面点较多的情况。首先，在点云空洞提取的基础上估计水面高程。然后，基于点云高程沿三角网进行边界扩张，直至边界不再变化。在扩张中，动态更新对水面高程的估计。最后，修剪边界。本文采用江苏常州小黄山地区的点云进行了实验，成功提取到所有主要静态水体的边界，且与数字正射影像(digital orthophoto map，DOM)中的水体边界套合良好。

1 静态水体水面点高程分析

静态水体多为湖泊、池塘，或流动缓慢的河流(不考虑潮汐河口或水体与海水联通的情况)，短时间内水位变化小。对一片联通水体来说，液态水作为一种流体，在无其他外力作用时，会在重力的作用下流动，直至水体表面没有重力势能差。对中小面积的静态水体，如果没有明显局部重力异常，可以简单认为，除出入水口外，此时水面的椭球体高程基本一致。

但是，波浪、水生植物、水质等因素可能造成静态水体水面高程不一致。波浪产生的水面点间的高程差一般在厘米级甚至分米级。在水生植物中，挺水植物可以形成明显高出水面的点云，而浮藻、浮叶植物、漂浮植物等产生的点高程十分接近水面高程。如果沉水植物离水面很近且水质不佳，则有可能水体不能完全吸收激光脉冲，沉水植物产生回波点，所以，点云中可能存在低于水面的点。不同类型的水生植物形成的点云差异较大，分析困难。但考虑到浮藻、浮叶植物、漂浮植物等的点高程接近水面高程，沉水植物的回波点少且容易被排除，而茂密的挺水植物经常遮盖水面，使人难以判断水体边界，所以在估计水面高程范围时，本文不考虑水生植物的影响。

假设LiDAR测量和波浪波动相互独立，则可以认为，水面高程的方差Vwater由两部分构成：LiDAR测量相对误差的方差Vmeas和波浪造成的方差Vwave，见式(1)。

Vwater=Vmeas+Vwave

(1)

2 方法

本文提出了一种沿Delaunay三角网迭代扩张，动态估计水面高程的静态水体边界提取方法(图1)。考虑到水体内的岛屿、设施等会使水体内部存在空洞，将单个静态水体抽象为OpenGIS简单特征规范[12]中的多边形(Polygon)，允许空洞存在。

图1 本文方法流程图

2.1 点云去噪

为避免高程异常点影响水面高程估计，对点云进行直通滤波，只保留高程位于合理范围内的点。

2.2 初始边界计算

在计算水体初始边界时，本文对文献[13]的方法进行了修改，抛弃点云外边界，只留下点云空洞边界。设边长阈值为Lthre，面积阈值为Athre。在得到点云空洞边界后，构建带拓扑关系的水体初始边界，规定外边界点逆时针排列，空洞边界点顺时针排列。

此处定义：设边界上有P0、P1、…、PN共N+1个顶点，P0=PN，则称PiPi+1(0≤i≤N-1)为边界段；如果PiPi+1是由ΔPiPkPi+1搜索到，则ΔPiPkPi+1的邻接三角形中，Pk对着的那个三角形为搜索方向。

2.3 水面高程估计

假设机载LiDAR对同一个高程的观测值服从正态分布。通过统计测区内局部平坦地面的点云高程，可以大致估计Vmeas和对应的标准差σmeas。Vwave难以实测，但与水面波动幅度正相关。根据式(1)，Vwater≥Vmeas，所以σwater≥σmeas，其中σwater是水面高程的标准差。

本文认为水体边界上含有水面点，可以由外边界和空洞边界估计水面高程范围。和落在水体周围的植被、建筑物、地面上的点相比，水面点高程相对偏低。极少量来自沉水植物的点，则明显低于水面点的正常高程范围。为确定水面高程的可能范围，采用聚集层次聚类[14]对边界点高程进行分类。当类间最小距离均大于ds时，聚类终止。对每个高程类，记录高程最小值、最大值、高程数和点数。

在高程分类后，对所有高程类按照最小值从小到大的顺序进行遍历。最先满足高程数大于等于Nzthre的高程类被认为含有水面点高程。如果没有高程类满足条件，可适当增大ds，进行二次分类。如果二次分类还没有满足条件的高程类，则认为该边界不是有效的水体边界，抛弃该边界。设满足条件的高程类的最小值为Zmin，最大值为Zmax。在[Zmin，Zmax]内，值越小，是水面高程测量值的可能性越大，所以设置水面高程范围下限Zwmin=Zmin。为确定水面高程范围上限Zwmax，采用如下方法。

步骤1：对[Zmin，Zmax]内的每一个高程值Zt(受限于存储精度，Zt数量有限)，计算高程在[Zmin，Zt]内的边界点的平均高程Zmean和高程标准差σt。

步骤2：取边界点数大于5且σt≥σwater的最小Zt为Zwmax。如果没有Zt满足条件，则取Zwmax=Zmax。设Zwmax对应的Zmean为Zwmean。

2.4 水体空洞与建筑物空洞的区分

在机载LiDAR点云中，除了水体造成的空洞外，还经常存在建筑物造成的空洞(图2)。后者往往一侧为屋顶，另一侧为更低的屋顶或地面。上述特征使得建筑物空洞的边界点高程分布直方图呈现明显的双峰特点。

图2 建筑物空洞示例

使用如下方法区分建筑物空洞和水体空洞。

步骤1：采用2.3节的方法对空洞边界点进行高程分类。

步骤2：寻找高程数最多的2个高程类。

步骤3：如果2个高程类的点数均超过空洞边界点总数的1/3，且较低高程类的最大值低于较高高程类的最小值2.5 m以上，则认为该空洞是建筑物空洞，否则是水体空洞。

2.5 基于高程的边界扩张

在得到水体初始边界后，对每个水体沿外边界和空洞边界进行基于高程的边界扩张。首先定义，如果某三角形三个顶点高程均在[Zwmin，Zwmax]内，则该三角形为水平三角形，否则为非水平三角形。在单次扩张中，沿着搜索方向试图扩张一个三角形的距离，将水平三角形扩张到边界内部。设当前对边界段PiPi+1进行扩张，方法如下(参考图3)。

步骤1：删除已扩张边界的最后一个顶点。

步骤2：计算顶点含Pi的三角形的集合SΔ，以Pi为中心逆时针排列。

步骤3：筛选SΔ中位于PiPi+1和已扩张边界的最后一段PlPi之间，而且包含搜索方向的三角形集合St，以Pi为中心顺时针排列。

注：白箭头表示搜索方向。在(b)～(i)中，黑箭头表示扩张后的有序边界；白三角形为水平三角形；灰三角形为非水平三角形。图3 基于高程的边界扩张示意图

图3(b)至图3(i)显示了在图3(a)St中的3个三角形的8种组合(是否水平三角形)对应的扩张边界。对所有边界段扩张一遍，并排除无效边界，称为一轮扩张。对边界进行多轮扩张，至边界不再变化，则基于高程的边界扩张终止。在每轮扩张后，都需要动态更新对水面高程的估计。

2.6 交叉多边形处理

在边界扩张的过程中，扩张后的多个边界之间可能发生交叉。例如，同一片水体被浅滩分为多个部分，每个部分的边界扩张到越过浅滩导致边界间发生交叉。单个边界可能发生自交叉，而自交叉多边形不符合OpenGIS简单特征规范。为此，需要对交叉多边形(包括多个多边形间交叉和单个多边形自交叉)求并集。图4为示意图。图中，粗闭合实线为一个凹多边形的边界(粗花箭头表示搜索方向，实粗箭头表示边界点顺序，虚粗箭头表示边界点逆序)；细闭合实线为一个带洞多边形的边界(细花箭头表示搜索方向，实细箭头表示边界点顺序，虚细箭头表示边界点逆序)；1～6是候选多边形编号。首先，对边界采用文献[15]的方法进行分割。

图4 交叉多边形处理示意图

步骤1：对交叉边界按照交点分割为正序边界段。

步骤2：对正序边界段生成对应的逆序边界段。

步骤3：使用最小角准则连接所有边界段(包括正序和逆序边界段)，生成的多边形称为候选多边形。

由图4可以发现，1的边界是并集外边界，搜索方向朝外，逆时针排列；4和7的边界是并集空洞边界，搜索方向朝内，顺时针排列。本文规定：边界点逆时针排列且搜索方向朝外的候选多边形的边界为并集外边界，边界点顺时针排列且搜索方向朝内的候选多边形的边界为并集的空洞边界。此处只保留这两种候选多边形的边界。

两帮劈裂示意类似于图1a与图1b，在巷道开挖前，两帮内有原生节理裂隙及采动产生的裂隙，开挖后在强卸载作用下节理裂隙很快受剪切、拉伸破坏并贯通，发生剪胀变形，如果受到动载冲击波冲击，则这一过程加剧。贯通后的节理裂隙继煤体深部延伸，剪胀变形加剧，直发生大变形甚至失稳。

2.7 孤立点边界移除

由于水面高程估计不一定非常精确，所以可能有少量水面点的高程不在[Zwmin，Zwmax]内。这导致在基于高程的边界扩张后，可能出现一些孤立点造成的边界(图5)。对这些孤立点产生的边界，需要将它们排除。这些边界的顶点相对较少，边界段搜索方向朝内，且边界内部由环绕孤立点的三角形构成。基于以上特征，使用如下方法进行处理。

注：白点为高程在[Zwmin，Zwmax]内的点；黑点为高程不在[Zwmin，Zwmax]内的点；粗黑线为基于高程扩张后的边界。图5 孤立点边界示意图

步骤1：如果某边界的边界段数目小于等于Eiso，且搜索方向朝内，则遍历边界段，进入步骤2。

步骤2：对边界段PiPi+1，设代表其搜索方向的相邻三角形为Δneib，查找在Δneib中除Pi和Pi+1外的另一个顶点Pv。如果Pv不在边界上，则将Pv放入集合SV。

步骤3：遍历完后，如果SV中不重复顶点的数目小于等于Niso，则认为边界由孤立点造成，抛弃该边界，否则保留该边界。

2.8 边界修剪

上述步骤得到的水体边界可能存在重合边界点(非首尾点重合)造成的突出部分(简称突部)和内陷部分(简称陷部)。为了让边界平滑，最好从水体主边界上去除细碎的突部和陷部。在图6中，水体0是外边界包围部分1减去陷部对应的空洞2，而1由水体主体3和突部4构成。简单地通过面积阈值即可判断是否该保留陷部和突部。

图6 边界修剪示意图

方法如下。

步骤1：对待修剪边界计算外边界和空洞边界。

步骤2：如果空洞面积小于Athre，则抛弃空洞。

步骤3：对外边界利用最小角准则进行连接，得到1个或多个不重叠的多边形，称为子多边形。

步骤4：如果子多边形面积小于Athre，则抛弃子多边形，否则保留。

步骤5：对剩余空洞和子多边形重构拓扑关系。

3 实验

本文利用2019年9月由RIEGL VQ-1560i机载LiDAR采集的江苏省常州市小黄山地区的数据进行实验。小黄山的地类以森林、农田、居民区、道路和水体为主。在测区内分布着200多个水体，绝大部分为静态水体，包括天然湖、鱼塘、沼泽等，最大面积约62 000 m2，最小面积约120 m2。数据采集时设计地面点密度为16点/m2。除点云外，还同时获取了可见光影像，生成了0.1 m分辨率的DOM。在DOM辅助下，选取了220份水体数据，每份数据包括一个或多个邻近水体。采用C++编程实现了本文方法。取Lthre=1.0 m，Athre=1.0 m2。

3.1 LiDAR测量误差分析

在DOM辅助下，选取了5块位于局部近似平面的点云，去噪，计算高程标准差作为点云测量的高程相对误差σmeas，结果见表1。可以看出，σmeas分布在[0.017，0.025]内。因此，σwater需要大于等于0.025 m。

表1 局部近似平面的高程标准差

3.2 水面高程范围估计

图7是典型的水体边界点高程分布直方图。在11.47、11.52和11.54 m处存在3个离群低点。而从11.56 m开始到11.98 m，每个厘米值上均有点分布。设ds=0.01 m，则高程在[11.47，11.98]内的点共分为4个高程类。取Nzthre=6，则Zmin=11.56 m，Zmax=11.98 m。

图7 典型的水体边界点高程分布直方图

图8展示了Zwmean和σwater对Zwmax的依赖关系。可以看出，Zwmean和σwater是Zwmax的单调递增函数。σwater越大，对应的Zwmax越大。

图8 Zwmean和σwater随Zwmax变化的折线图

3.3 水面高程标准差对边界提取的影响

σwater取值影响水面高程范围估计，进而影响最终结果。为此，实验了0.026、0.028、0.030、0.032和0.034 m共5个σwater值用于水体边界提取的效果，取Eiso=12，Niso=2。基于不同σwater得到的水体边界一般存在差异。对结果准确度进行定量评估非常困难。但是，通过与测区DOM叠加显示，可以目视判断由点云提取的水体边界和影像中边界的套合程度，作为评估结果的依据。

结果表明，对所有220个数据来说，这5个σwater值均100%提取到主要水体。在水体内没有很多点的情况下，即使σwater较小也能取得良好的效果，使用较大的值并不会明显改进效果。

图9的鱼塘中，σwater取0.026 m和0.034 m得到的边界基本一致，只有细微差别。

但是，对于水体内存在大量点的情况，σwater可能明显影响边界提取的效果。

图10中的4个完整鱼塘，除了中间的小鱼塘外，其他三个大鱼塘在中南部均存在大量点，其中极少量点位于增氧机，其余点位于水面。这可能是在采集数据时鱼塘南部刮南风造成的，未影响到鱼塘北部。在σwater=0.026 m时，三个大鱼塘内部都存在大量空洞。随着σwater的增大，空洞越来越少。到σwater=0.034 m时，除了鱼塘内的增氧机外，已经没有任何水面点造成的空洞，而鱼塘外边界、增氧机边界与DOM中的边界套合良好。

注：在(b)、(c)内，黑色实线为提取到的边界。图9 水体内点较少情况时的边界提取效果示例

注：在(b)至(f)内，黑色实线为提取到的边界。图10 水体内点较多情况时的边界提取效果示例

4 结束语

水面高程标准差σwater的取值对水体内点少的情况影响小，但对水体内点较多的情况影响大，较大的σwater倾向于给出较好的效果。对后者来说，增大σwater会增大Zwmax，使更多波浪点的高程小于等于Zwmax，从而在基于高程的边界扩张中被扩张到。对一些动态水体，比如喷泉，采用较大的σwater也可以取得良好效果。但是，对于建筑物和水体相邻的情况，如果建筑物空洞与水体空洞相连，则提取到的水体边界可能不准。本文方法不适用于没有足够大面积的空洞存在，以及挺水植物较多的水体。

本文方法提取的水体边界可以辅助测绘内业中DEM水体断裂线的生产，一般经少量人工编辑即可获得较常规生产手段更真实准确的水体边界，大大减少人工工作量。

水体边界提取是空间不确定性研究中的典型问题[16]。由于水体自身特点和观测手段的限制，难以精确获取某时刻水体的边界。如何定量评估水体边界的准确度需要进一步研究。