在机载LiDAR数据中, 纹理可定义为因局部区域高程变化而产生的对比度、均匀性等物理特征, 即高程纹理^[6]。由图2可以看出, 建筑物与邻近点的高差比较固定, 没有明显的高程突变(图2的P1, P2点), 可以近似看作连续的水平平面; 植被间的高差较大, 且不规则(图2的P3, P4, P5点), 其与邻近点之间的斜率很大, 很难看作近似的水平面。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 LiDAR点云空间分布Fig.2 Spatial characteristics of LiDAR point clouds

为了提高建筑物点云提取的精确性, 本文在提取建筑物点云的过程中加入邻近点的斜率信息, 通过设定一定的参数来拟合这种近似平面(图3)^[7]。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 基于高程纹理提取建筑物点示意图Fig.3 Detecting building points based on height texture

图3上P点为建筑物中心点; P₁为相邻点; P10为P₁到中心点P所在平面的投影点; S为P₁P之间的空间距离; d为P₁P之间的平面距离; h为P₁P之间的高差; θ 为P₁点与P点所在平面的夹角。设定高差阈值H, 平面距离阈值R和角度阈值⌀, 判定P₁与P是否在同一个平面上需要满足3个条件, 即d< R且|h|< H, θ < ⌀。

一般高差阈值H和平面距离阈值R分别设为平均高差和平均平面距离的1~2倍; 倾角阈值⌀设为10° ~20° 。拟合倾斜平面所需要的参数阈值一般大于平面的参数阈值。过小的阈值会增大平面点的遗漏误差, 过大的阈值将增大平面的侵蚀作用。为了避免某些少量满足条件的点被归为平面, 需设定阈值N, 只有满足条件的P点所有邻近点的数目大于N, 才能归为平面, 认为P点为建筑物点; 反之, 为非平面, 即P点为其他地物点。

建筑物边缘特征是构成建筑物模型的重要信息, 能够描述建筑物结构和几何特征, 大多数建筑物的特征线均为直线, 检测到建筑物的各个直线特征就意味着检测到了建筑物的主要结构信息。因此, 建筑物直线特征的检测是建筑物提取的重要工作。一般常用的边缘检测算子有梯度算子、Laplace算子、Sobel算子、Canny算子以及方向算子等非线性算子, 此外还有曲面拟合法等等。其中, Canny算子具有方向性, 能更好地应用于边缘强度估计, 能产生梯度方向和强度2个信息^[8]。因此, 本文利用Canny算子检测建筑物脚点深度影像中的建筑物边缘。

本文采用Ruijin ^[9]提出的方位角聚类比较法规则化建筑物轮廓。该方法结合了聚类和调整2个过程, 聚类方法类似于K-means方法。其算法描述如下:

1)计算建筑物线段各自的方位角, 根据方位角将线段分为2类。求这2类方位角的平均值a和b, 比较候选线段与这2个平均值的差别, 差异小的为该线段所属的类别。经处理后, 这2类线段间应为互相垂直关系。

2)在每一类中计算方位角的权重均值, 因为较长线段比较短线段方位角精度高, 所以权重取决于线段的长度, 即

α _azim=∑ l_i· αazimi/∑ l_i(1)

式中: l_i为其中一类中第i条线段的长度(i=1, 2, 3, …, n); αazimi为第i条线段的方位角; α _azim为其中一类的方位角。

3)使用Gauss-Markov模型调整2类权重, 使得2类线段严格垂直。这里的权重为这2条线段的总长度。

4)将每一条线段绕该线段中点旋转至该类方位角。到此为止, 建筑物的轮廓线段之间不是垂直, 就是平行。

5)设定平行线段间距离阈值d(可取1 m), 当2条线段间距离小于该值, 则融成一条新线段。通过计算各线段中点的权重均值获取新线段的中心点。比如, 当2条线段融合成新的线段时, x坐标求取公式为

x= (l1x1+l2x2)l1+l2, (2)

式中: l₁和l₂分别为2条线段的长度; x₁和x₂为2条线段的中心点坐标。

该算法的优势在于无需较好的建筑物初始轮廓, 就能获取规则的建筑物边界。也就是说, 不需要非常准确的边缘检测算子即可提取出建筑物轮廓。

对于约束D-三角网剖分算法, 根据约束边嵌入时机的不同可分为2类: 第1类是指在构网的同时考虑约束边的影响, 直接构建约束D-三角网; 第2类是指首先不考虑约束边的影响, 构建数据集的非约束D-三角网, 在已构好的三角网中强行嵌入约束边以构建约束D-三角网, 即2步法^[10]。本文采用2步法进行构网, 嵌入算法的步骤如下:

1)按照一定规则对数据区域进行格网划分, 基于逐点内插算法对离散地面点集生成D-三角网。

2)确定线段的影响区域。如图4所示, 首先确定线段AB的端点A所在的任一△ A, 然后由△ A来确定线段方向上的首三角形(记为△ S); 在△ A(如△ APB)中, 确定点A的对应边b(如 PB)与有向线段AB相交否, 若相交, 则首三角形为△ A(△ APB=△ S); 如果判断边c(如DE)在有向线段AB的右侧, 则以端点A为圆心逆时针方向搜索与△ A相邻的下一个三角形; 若为左侧, 则以顺时针方向搜索。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 搜索线段影响区域Fig.4 Detecting affecting area of the line

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 影响区域的三角剖分Fig.5 Delaunay triangulation by incremental insertion algorithm for affecting area

3)影响区域的三角剖分。如图5所示, 以有向线段AB作为扩展边, 在扩展边右侧影响区域(图5中阴影部分)的点集中取一点C, 使得C点与扩展边的两端点的连线组成的夹角最大(即最大角准则), 生成新的△ ABC。同时, 在有向线段AB左侧构网时, 则需要将有向线段AB改为BA。将多边形的每条线段都重复上述过程, 最终实现多边形入网。

4)多边形内部清空处理。从多边形内部出发, 根据拓扑关系, 向八方向辐射, 搜索位于多边形内部的三角网并将其移除。图6(a)为待移除位于内部多边形的三角网; 图6(b)为重新构造的D-三角网。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 嵌入约束多边形前(左)后(右)D-三角网的变化示意图Fig.6 Comparison of constrained delaunay triangulation(left) and delaunay triangulation(right)

实验区位于广西壮族自治区柳州城区, 其正射影像如图7所示。区内有多种地物, 如植被、道路、汽车及大小高度形状各异的建筑物等。实验数据由ALS50-II系统获得, LiDAR点云密度为8点/m²左右, 数据的垂直精度优于15 cm, 水平精度优于0.5 m, 总点数为1 976 040。图8为实验区激光点云高程设色图。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 正射影像Fig.7 Orthophoto

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 高程设色图Fig.8 Display at height

首先, 对点云进行预处理。利用TerraScan滤波得到地面点数据, 实验区属于丘陵地区, 角度阈值设为6.0° ; 根据测区的实际建筑物面积, 将建筑物最大边长设为160 m; 构建三角形过程中的高差阈值设为1.4 m。当所加点构成的三角形每条边短于5 m时, 阻止向三角形内部加点。滤波后获取的地面点如图9所示。可以看出, 对于植被和建筑物等高程突变比较明显的地物来说, TerraScan滤波方法具有比较好的效果。

其次, 进行建筑物点的检测。首先获取去除地面点的nDSM(normalized DSM)(图10), 利用3 m的高度阈值滤掉低矮灌木丛, 为下一步的建筑物平面拟合减少误差; 然后基于高程纹理提取建筑物点, 令阈值R=4 m, H=1 m, ø=15° , N=40。建筑物点的提取结果如图11所示。

然后, 提取建筑物矢量轮廓线。首先, 基于Canny边缘检测算法提取由建筑物点云所生成的nDSM深度图像中建筑物边缘(图12); 然后, 利用方位角聚类比较法规则化建筑物轮廓(图13); 最后, 在规则化的建筑物二维轮廓线上任取一点, 由离散三维点云内插出其高程, 将高程值赋与此建筑物。应用该方法得到的建筑物轮廓线与本区域的正射影像叠合如图14所示。可以看出, 提取的建筑物边缘基本与正射影像上的建筑物边缘重叠。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 地面点Fig.9 Ground points

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 滤除地面点的nDSMFig.10 nDSM after filting ground points

图11

Fig.11

	Figure Option ViewDownloadNew Window
	图11 建筑物点云Fig.11 Building points

图12

Fig.12

	Figure Option ViewDownloadNew Window
	图12 Canny算子提取的建筑物边缘Fig.12 Building edges detected by Canny

图13

Fig.13

	Figure Option ViewDownloadNew Window
	图13 建筑物边缘规则化Fig.13 Building edges normalized

图14

Fig.14

	Figure Option ViewDownloadNew Window
	图14 建筑物轮廓与DOM叠合Fig.14 Combination of building contour and DOM

最后, 利用滤波得到的离散三维地面点与建筑物矢量轮廓线构建TIN格网, 这里格网间距设为1 m, 得到DSM的深度图像(图15)。

图15

Fig.15

	Figure Option ViewDownloadNew Window
	图15 DSM深度图像Fig.15 Depth image of DSM

如图14所示, 实验区域内共有71栋建筑物, 其中有66栋被检测出来, 对于建筑物顶部为平面的规则建筑物检测效果比较明显, 而对于不规则和尖顶房有漏检的现象。将提取的结果与通过TerraScan软件(邻近点高差阈值设为0.2 m)自动提取的建筑物结果相比较(图16, 17), 可以看出, 本文方法对于建筑物的提取更为精确, 建筑物边缘有很好的保留, 几乎没有建筑物边缘点被误分为高植被点的现象。

图16

Fig.16

	Figure Option ViewDownloadNew Window
	图16 TerraScan提取的建筑物点和高植被点Fig.16 Building points and high vegetation points detected by TerraScan

图17

Fig.17

	Figure Option ViewDownloadNew Window
	图17 本文方法提取的建筑物点和高植被点Fig.17 Building points and high vegetation points detected by the paper method

为了定量描述建筑物的分类误差, 本文以TerraScan的半自动、半手工分类结果为参考数据, 计算本文方法的分类误差(表1)。实验结果表明, 本文方法的分类误差小于10%, 说明方法比较精确有效。

表1

Tab.1

表1(Tab.1)

表1 分类误差分析 Tab.1 Error analysis of classification 参考数据/个 本文结果/个 误分点/个 误分点比例/% 792 688 728 520 641 68 8.09

表1 分类误差分析 Tab.1 Error analysis of classification

比较本文方法构建的DSM(图18)及采用逐点内插得到的规则格网(格网间距也为1 m)所构建的DSM(图19)建筑物细节, 可以得出, 本文构建的DSM建筑物边缘更准确和规则。

图18

Fig.18

	Figure Option ViewDownloadNew Window
	图18 本文方法构建的DSM建筑物边缘Fig.18 Building edge of DSM extracting by proposed method

图19

Fig.19

	Figure Option ViewDownloadNew Window
	图19 规则格网构建的DSM建筑物边缘Fig.19 Building edge of DSM based on regular grid

建筑物矢量轮廓线的高程信息是由二维边缘线上任意一点内插得到的, 所以构建的DSM存在一定的误差。本文利用抽稀建筑物激光点云内插得到的DSM高程误差D_Z如表2所示。高程平均误差D_Z平均=0.018 5 m, 均方根误差RMSE=0.069 5 m。表明利用本文方法所构建的DSM建筑物边界信息比较精确, 高程精度比较高。

表2

Tab.2

表2(Tab.2)

表2 DSM高程误差 Tab.2 Evaluation of DSM accuracy X/m Y/m Z/m DZ/m 334 635.44 2 69 0427.09 94.26 0.051 334 607.60 2 69 0427.13 91.58 0.011 334 615.75 2 69 0427.17 91.62 0.000 334 634.33 2 69 0427.22 94.46 0.012 334 608.95 2 69 0427.29 93.67 0.033 334 626.53 2 69 0427.34 91.59 0.003 ︙ ︙ ︙ ︙ 334 637.63 2 69 1067.10 91.74 0.015 334 637.61 2 69 1138.59 90.84 0.007 DZ平均=0.018 5 RMSE=0.069 5

表2 DSM高程误差 Tab.2 Evaluation of DSM accuracy