基于城墙断面从LiDAR点云数据自动提取古城墙数据的方法

摘要

本发明涉及一种基于城墙断面从LiDAR点云数据自动提取古城墙数据的方法，包括以下步骤：获取城墙待选区，矢量化得到顶部二维中心线，并插值得到完整的城墙顶部三维中心线；对城墙顶部三维中心线密集采样得到探测基点，计算每个探测基点断面上的点云；构建城墙断面格网；并基于城墙断面“横向探测—纵向探测”城墙顶部、垛口和侧面点云数据。本发明提出利用LiDAR点云数据提取古城墙的数据，解决了利用LiDAR技术进行数字存档的关键技术，有效的解决了传统测量耗时耗力，且无法获取城墙侧面细节信息的不足。本发明提出的“纵向探测‑横向探测”的城墙侧面点精确提取方法，可以有效的实现城墙侧面细节信息和顶部信息的自动获取。

说明书

技术领域

本发明涉及一种提取城墙点云数据的方法，特别是一种基于城墙断面从LiDAR点云数据提取城墙点云数据的方法。

背景技术

古城墙是中国文化遗产的重要组成部分，有着重要历史和文化价值。但由于人为的破坏、战争的摧毁、自然的风化、洪水的毁坏等，古城墙正遭受严重的破坏。对古城墙进行数字存档是重要的保护手段，这不仅有利于了解古城墙的现状，而且对于古城墙的复建和保护有着重要的指导意义，也对于古城墙的展示和宣传有着重要作用。

城墙构造的目的是用于城市的抵御，通常依据城池的布局而建，在空间上蜿蜒几到几十公里；为了抵御敌人的攻击，往往建为城墙基部较厚、顶部较窄，呈现“上窄下宽”(如图1所示城墙断面)，且侧面呈现非线性变化的结构特点。目前获取城墙信息，大多采用人工测量的方式，十分耗时耗力。由于测量技术的制约，人们对城墙的认识和了解一直不够全面，历次城墙“现状”调查中，各种数据亦多有出入。现有的测量方式在古城墙的数字存档中有很多不足，具体可以归纳为以下两点：

1)测量中采用传统测量和影像量测相结合的方式，传统测量只能测得离散的点、线信息，无法对城墙“上窄下宽”不规则变化的侧面细节进行面覆盖的数据采集，影像量测由于卫星影像的分辨率和二维特性的制约，也无法实现侧面的信息采集；

2)采用分段量测城墙的高度，无法完整地获取城墙整体的高度，不利于掌握城墙的高度变化，不利于城墙信息的及时更新和数字存档。

激光雷达技术是近年来快速发展的一种新型测量技术，能够大范围、高时效、无接触式地获取空间信息，已经成为辅助文化遗产保护的重要手段。通过机载激光雷达对城墙进行扫描，得到的点云数据包括地面点、城墙点、一般建筑物点、植被等。从上述原始点云中自动提取城墙点云数据，是古城墙数字存档中的关键技术。利用机载激光雷达数据自动提取城墙点云存在以下难点：

1)“上窄下宽”且不规则变化的侧面造成点云精确提取困难。一般建筑物的侧面与地面垂直，从机载激光雷达点云数据中提取建筑物，常用的方法有边缘检测、区域生长、数学形态学等，这些方法的实质是依赖建筑物的高程突变，得到建筑物的顶部点云或者顶部轮廓。而城墙点云提取包括顶部和侧面点云数据的提取，但城墙“上窄下宽”的结构使得其侧面点云与地面点云在交界处，高程突变不明显，而且侧面呈现非线性的不规则变化，因此难以用现有的算法提取城墙点云。

2)如图1所示，城墙点云数据存在缺失且有其他干扰。点云数据的缺失主要由遮挡引起，遮挡的类型有三类：一是由于航线引起的城墙自身遮挡，二是城墙附近植被的遮挡，三是城墙附近棚户对城墙侧面的遮挡。遮挡造成城墙侧面数据不完整，而往往在缺失侧面数据处，又通常存在紧邻的植被或者棚户和房屋(如图1所示)，造成对城墙点云数据的干扰。数据的缺失和其他地物的干扰，进一步增加了城墙侧面提取的难度。

发明内容

本发明要解决的技术问题是：提出一种从LiDAR点云数据中提取古城墙数据的自动高效的方法，该方法可以克服原始数据中数据缺失和其他地物的干扰，能够完整呈现城墙侧面的非线性不规则变化，从而实现对古城墙点云数据的数字存取。

本发明为解决上述技术问题提出的技术方案是：一种基于城墙断面从LiDAR点云数据自动提取古城墙数据的方法，包括以下步骤：

第一步、提取城墙待选区

获取预定区域的机载LiDAR点云数据，从所述LiDAR点云数据中分离出地面点云数据，得到该预定区域的DEM数据和DSM数据，通过栅格计算从DSM数据中减去DEM数据得到归一化的DSM数据，即得到nDSM数据；从所述nDSM数据中滤除城墙附近的建筑物和植被，得到多个连通区域，其中面积最大的连通区域即为城墙待选区；

第二步、从所述城墙待选区中找出城墙顶部

构建三维体元坐标系，得到所述城墙待选区的三维体元格网，记录所有体元内部所含的LiDAR数据点，剔除没有LiDAR数据点的空体元后计算剩余每个体元的粗糙度；遍历所述城墙待选区所有的体元，若以每个体元为中心的3×3×3空间邻域中的相邻体元的粗糙度差值小于预设的第一粗糙度差值阈值，则该中心体元与相邻体元内的LiDAR数据点属于同一类别；计算各类别的体元中所有LiDAR数据点的平均高程值，其中平均高程值大于预设的高程阈值且体元个数最多的一类体元即为城墙顶部；

每个体元的粗糙度的计算方法如下：假设一个体元中共有n个LiDAR数据点，h_i为该体元中第i个LiDAR数据点的高程值，其中1≤i≤n；则该体元的粗糙度式中

第三步、提取城墙顶部的三维中心线

将城墙顶部的点云数据转换为面状要素，然后对所述面状要素进行细化和跟踪并矢量化后得到城墙顶部的二维中心线；对所述二维中心线进行等间距重采样得到多个采样点，按delaunay三角剖分算法找到每个采样点周围构成三角形的属于城墙顶部的三个LiDAR数据点，然后对上述三个LiDAR数据点的高程值进行线性内插得到每个采样点的高程值；将每个采样 点赋予高程值后再矢量化从而得到城墙顶部的三维中心线；

第四步、构建城墙断面空间格网

将所述城墙待选区的LiDAR点云数据和所述三维中心线均投影至一个二维平面，沿所述三维中心线投影至该二维平面得到的城墙中心线按小于平均点间距的距离进行等间距重采样得到m个探测基点，依次是O₁，O₂，…，O_m，m≥2，所述平均点间距是机载LiDAR点云数据在扫描测量时的采样间隔，是数据本身的一个属性，可以通过查看原始LiDAR数据获得；对每个探测基点O_I，1≤I≤m，计算与探测基点O_I邻近的探测基点O_I-1构成的直线L1，该直线L1为探测基点O_I的探测基准线，与探测基准线L1垂直且交于O_I点的垂线L2即为探测基点O_I的断面；遍历投影至所述二维平面的城墙待选区所有LiDAR点云数据，每个到垂线L2的距离小于为断面厚度length的LiDAR数据点均属于所述探测基点O_I的断面；构建探测基点O_I的断面空间格网，令城墙断面体元的宽度width等于平均点间距，城墙断面体元的高度height为所述城墙顶部的LiDAR数据点的最大高度和最小高度之差；计算探测基点O_I的断面上的每个LiDAR数据点到直线L1的平面距离d1，并根据该LiDAR数据点在直线L1的不同侧赋予flag值为1或-1，从而得到该LiDAR数据点的行、列、层号(e,j,k)＝(I,Int(flag*d1/width),Int((z-Z_I)/height)+Number)，式中e是该LiDAR数据点在城墙断面体元格网中的行号，j是该LiDAR数据点在城墙断面体元格网的中的列号，k是该LiDAR数据点在城墙断面体元格网中的层号，number是用于保证k≥0的常数，z是该LiDAR数据点的高程值，Z_I是探测基点O_I的高程值；

第五步、基于城墙断面格网提取城墙顶部、垛口和城墙侧面

对每个探测基点O_I进行横向探测，以探测基准线L1为基准在城墙断面格网内向两侧探测，若探测的体元内有LiDAR数据点，则该体元为顶部体元，计算该体元的粗糙度，然后计算此体元与后一体元的粗糙度差值，若该粗糙度差值小于预设的第一粗糙度差值阈值，则后一体元也为顶部体元；若该粗糙度差值大于预设的第二粗糙度差值阈值，则后一体元为垛口体元，并向上探测下一体元，若下一体元内有LiDAR数据点，则该下一单元也为垛口体元；当探测到垛口体元后向外探测一个体元，若该体元内有LiDAR数据点，则该体元为顶部体元；若在同一个方向连续探测到两个或两个以上没有LiDAR数据点的空体元时，则到达顶部边缘，停止该方向扫描；

对每个探测基点O_I进行纵向探测，以探测基点O_I所在层为基准往下移动一层，以纵向探测基准线L1为基准向两侧扫描探测，直至探测次数超过预定次数或者扫描到体元内有LiDAR点数据，则停止该侧扫描，内有LiDAR点数据的体元为城墙侧面；继续往下移动一层，重复向两侧扫描探测，直到探测基准线L1两侧的体元内均没有LiDAR数据点，则停止纵向探测。

本发明中DEM为数字高程模型(Digital Elevation Model)，DSM为数字表面模型(Digital Surface Model)，两者均为常用的地理信息数据。本发明第一步中从所述LiDAR点云数据中分离出地面点云数据后，提取得到预定区域的DEM数据和DSM数据为现有技术，可以采用多种现有方法实现。第三步中delaunay三角剖分算法也为现有技术。

优选的，第一步中采用不规则三角网加密滤波算法从所述LiDAR点云数据中分离出地面点云数据。

优选的，第一步中所述nDSM数据中滤除城墙附近的建筑物和植被的方法是：对所述nDSM数据的高程值进行二值化操作，滤除高程值低于所述预设的高程阈值的nDSM数据即滤除了城墙附近的植被和建筑。

由于城墙侧面高度渐变，与地面分界不明显，为了保证城墙点云数据的完整性，对第一步中得到的城墙待选区进行缓冲操作，即对所述城墙待选区向外扩充预定距离。

本发明带来的有益效果如下：

1)本发明利用LiDAR点云数据提取古城墙的点云数据，可以方便对古城墙进行数字化存档，有效的解决了古城墙传统测量方式耗时耗力，且无法获取城墙侧面细节数据(信息)的不足。

2)本发明通过“纵向探测-横向探测”提取精确城墙侧面点云数据，可以有效的实现城墙侧面细节信息的自动获取，解决了传统测绘方式因城墙侧面不规则变化导致的难以获得城墙侧面细节信息的难题。

3)本文提出的“城墙待选区提取-中心线获得-顶部侧面点云提取”技术方案，可以有效从大区域机载LiDAR点云数据中，可以非常准确地分离出城墙侧面和顶部的完整信息，解决了利用机载LiDAR点云数据进行古城墙数字存档的关键难题。

4)本发明方法适应性强，经实践证明，使用本发明方法以南京城墙为例的古城墙点云数据进行提取，与人工提取结果相比，其提取正确率和完整率均达97％以上，完全可以实现自动化提取古城墙的点云数据。

附图说明

下面结合附图对本发明进一步说明。

图1是一般古城墙的断面示意图。

图2是本发明实施例待测实验区的LiDAR点云数据示意图。

图3是本发明实施例中试验区A提取的城墙数据示意图。

图4是图3中a部的放大示意图。

图5是图3中b部的放大示意图。

图6是图3中c部的放大示意图。

图7是本发明实施例中试验区B提取的城墙数据示意图。

图8是图7中d部的放大示意图。

图9是图7中e部的放大示意图。

图10是本发明实施例中试验区C提取的城墙数据示意图。

图11是图10中f部的放大示意图。

图12是图10中g部的放大示意图。

图13是本发明实施例中提取城墙顶部三维中心线并构建城墙断面格网的方法示意图。

图14是本发明实施例中对每个探测基点进行横向探测和纵向探测的方法示意图。

图15是本发明实施例的流程示意图。

具体实施方式

实施例

如图2所示，本实施例的实验数据是江苏省南京市城区内的机载LiDAR点云数据，其平均点间距为1m，高程精度为20cm，平面精度为50cm。

本实施例选择图2中城墙上的三段实验区域进行验证，参见图2中A、B、C所示区域，对LiDAR点云数据按高程进行渲染。第一段城墙实验区域A覆盖面积1300m*1000m，包含5589094个LiDAR数据点，城墙全长2062m，侧面数据缺失较多，城墙顶部边缘附着一定植被，另外侧面有城门结构。第二段城墙实验区域B，覆盖面积630m*630m，包含1674059个LiDAR数据点，城墙全长674m，城墙数据整体较为完整，城墙周围附着的植被和房屋较少。第三段城墙实验区域C，覆盖面积1500m*630m，包含1127355个LiDAR数据点，城墙全长1544m，侧面数据较为完整，侧面有城门结构，顶部和侧面边缘附着大量植被，且城墙顶部上有少量茂盛的植被。这三段墙实验区域几乎包含了城墙点云数据的所有特点，非常具有代表性。

本实施例利用LiDAR云数据对古城墙进行信息采集，并且对获取图2中所示的三段城墙实验区域A、B和C都做基于城墙断面的点云自动提取，方法基本相同。

本实施例以实验区域A为例进行说明基于城墙断面从LiDAR点云数据自动提取古城墙数据的方法，如图15所示，具体包括以下步骤：

第一步、提取城墙待选区

获取图2中所示实验区域A机载LiDAR点云数据，从所述LiDAR点云数据中分离出地面点云数据，得到该预定区域的DEM数据和DSM数据，通过栅格计算从DSM数据中减去DEM数据得到归一化的DSM数据，即得到nDSM数据；从所述nDSM数据中滤除城墙附近的建筑物和植被，得到多个连通区域，其中面积最大的连通区域即为城墙待选区。

本实施例对城墙实验区域A的LiDAR点云数据使用基于不规则三角网加密滤波算法实现 地面点云数据和非地面点云数据的分离，得到城墙实验区域A的地面点云数据。计算城墙实验区域ADSM和DEM，用DSM减去DEM，得到nDSM数据，也就是得到非地面物体。非地面物体包括城墙、一般建筑和植被等，有的建筑物之间，建筑物与植被之间可能存在连接。这里设置一个高程阈值进行二值化操作以滤除较低的植被和低矮的建筑，使得连接起来的建筑物被分开，得到一个二值化的建筑物标记影像。在标记影像中，将所有连通的区域标记为候选的建筑物区域并计算面积，筛选出面积最大的连通区域，该连通区域就是我们所需要的城墙待选区。

由于城墙侧面高度渐变，与地面分界不明显，为了保证城墙点云数据的完整性，对所述城墙待选区进行缓冲操作，即对所述城墙待选区向外扩充预定距离，缓冲操作可以很方便地通过现有软件Arcgis 10.1实现。

第二步、从所述城墙待选区中找出城墙顶部

以所述城墙待选区中的LiDAR点云数据的坐标最小值为原点并设定体元大小来确定，构建规则三维体元坐标系统，规则体元坐标系统中的体元是一个标准的立方体，得到所述城墙待选区的三维体元格网，记录所有体元内部所含的LiDAR数据点，剔除没有LiDAR数据点的空体元后计算剩余每个体元的粗糙度；遍历所述城墙待选区所有的体元，若以每个体元为中心的3×3×3空间邻域中的相邻体元的粗糙度差值小于预设的第一粗糙度差值阈值，则该中心体元与相邻体元内的LiDAR数据点属于同一类别；计算各类别的体元中所有LiDAR数据点的平均高程值，其中平均高程值大于预设的高程阈值且体元个数最多的一类体元即为城墙顶部。

每个体元的粗糙度的计算方法如下：假设一个体元中共有n个LiDAR数据点，h_i为该体元中第i个LiDAR数据点的高程值，其中1≤i≤n；则该体元的粗糙度式中

第三步、提取城墙顶部的三维中心线

将城墙顶部的点云数据转换为面状要素，然后对所述面状要素进行细化和跟踪并矢量化后得到城墙顶部的二维中心线；对所述二维中心线进行等间距重采样得到多个采样点，按delaunay三角剖分算法找到每个采样点周围构成三角形的属于城墙顶部的三个LiDAR数据点，然后对上述三个LiDAR数据点的高程值进行线性内插得到每个采样点的高程值；将每个采样点赋予高程值后再矢量化从而得到城墙顶部的三维中心线。

本实施例通过ArcGIS软件平台，将得到的城墙顶部LiDAR点云数据转换为面状要素，然后通过ArcScan的扩展模块，进行细化和跟踪等操作，完成城墙顶部二维中心线的自动提取。

第四步、构建城墙断面空间格网

如图13所示，将所述城墙待选区的LiDAR点云数据和所述三维中心线均投影至一个二维平面，沿所述三维中心线投影至该二维平面得到的城墙中心线按小于平均点间距的距离进行等间距重采样得到m个探测基点，依次是O₁，O₂，…，O_m，m≥2，所述平均点间距是机载LiDAR点云数据在扫描测量时的采样间隔，是数据本身的一个属性，可以通过查看原始LiDAR数据获得；对每个探测基点O_I，1≤I≤m，计算与探测基点O_I邻近的探测基点O_I-1构成的直线L1，该直线L1为探测基点O_I的探测基准线，与探测基准线L1垂直且交于O_I点的垂线L2即为探测基点O_I的断面；遍历投影至所述二维平面的城墙待选区所有LiDAR点云数据，每个到垂线L2的距离小于为断面厚度length的LiDAR数据点均属于所述探测基点O_I的断面；构建探测基点O_I的断面空间格网，令城墙断面体元的宽度width等于平均点间距，城墙断面体元的高度height为所述城墙顶部的LiDAR数据点的最大高度和最小高度之差；计算探测基点O_I的断面上的每个LiDAR数据点(如图13中p点)到直线L1的平面距离d1，并根据该LiDAR数据点在直线L1的不同侧赋予flag值为1或-1，从而得到该LiDAR数据点的行、列、层号(e,j,k)＝(I,Int(flag*d1/width),Int((z-Z_I)/height)+Number)，式中e是该LiDAR数据点在城墙断面体元格网中的行号，j是该LiDAR数据点在城墙断面体元格网的中的列号，k是该LiDAR数据点在城墙断面体元格网中的层号，number是用于保证k≥0的常数，z是该LiDAR数据点的高程值，Z_I是探测基点O_I的高程值；

第五步、基于城墙断面格网提取城墙顶部、垛口和城墙侧面

如图14所示，将城墙断面格网内的体元分为顶部单元和一般体元。对每个探测基点O_I进行横向探测，以探测基准线L1为基准在城墙断面格网内向两侧探测。以左侧为例若探测的体元内有LiDAR数据点，则该体元为顶部体元，计算该体元的粗糙度，然后计算此体元与后一体元的粗糙度差值，若该粗糙度差值小于预设的第一粗糙度差值阈值(第一粗糙度差值阈值，应略小于城墙顶部体元粗糙度和城墙顶部边缘垛口体元粗糙度差值)，则后一体元也为顶部体元；若该粗糙度差值大于预设的第二粗糙度差值阈值，则认为探测到达城墙顶部的垛口处，体元粗糙度出现突变，即后一体元为垛口体元，并向上探测下一体元，若下一体元内有LiDAR数据点，则该下一单元也为垛口体元，标记垛口体元为了获取完整的城墙点云数据，避免图中顶部边缘体元中数据的丢失；当探测到垛口体元后向外探测一个体元，若该体元内有LiDAR数据点，则该体元为顶部体元；若在同一个方向连续探测到两个或两个以上没有LiDAR数据点的空体元时，则认为扫描达到顶部边缘，此处垛口数据缺失，即到达顶部边缘，停止该方向扫描。右侧的扫描探测与左侧相同。

对每个探测基点O_I进行纵向探测，以探测基点O_I所在层为基准往下移动一层，以纵向探测基准线L1为基准向两侧扫描探测，直至探测次数超过预定次数或者扫描到体元内有LiDAR>

完成以上五个步骤后，就可以得到完整的实验区域A的城墙点云数据，包括顶面、侧面、垛口等，如图3-6所示。

采用同样的方法可以从实验区域B提取出如图7-9所示的城墙点云数据，以及从实验区域C提取出如图10-12所示的城墙点云数据。

以手动提取的城墙点云数据作为真实城墙点云数据，将本实施例的实验城墙点云数据点云和真实点云真实城墙点云数据进行叠合，获得两者的重叠点(正确提取的点云数据)、手动提取独有而实验结果缺失的点(遗漏点)以及实验结果独有而手动提取缺失的点(错误点)，点云提取的正确率和完整率计算公式如下：

其中，TP(True Positive)是指正确提取的城墙点云数据；FN(False Negative)是指遗漏的城墙点云数据；FP(False Positive)是指误分为城墙点的非城墙点，即错误提取的城墙点云数据。分别统计这三类点云数据的点数，按照上述公式计算正确率和完整率。

实验区域A的城墙点云数据提取结果的完整率和正确率为97.12％，98.23％，实验区域B的城墙点云数据提取结果的完整率和正确率为99.15％，98.97％，实验区域C的城墙点云数据提取结果的完整率和正确率为97.02％，98.28％。由上可知，采用本发明提取古城墙数据的方法提取结果的完整率和正确率均达到97％以上，这表明本发明利用LiDAR点云数据对古城墙进行信息采集并基于城墙断面自动提取城墙点云的方法可行，且有着较高的精度。

本发明不局限于上述实施例，本发明的上述各个实施例的技术方案彼此可以交叉组合形成新的技术方案，另外凡采用等同替换形成的技术方案，均落在本发明要求的保护范围内。