1 引言

三维城市模型的应用越来越广泛,涉及到生产生活等各个领域。该模型的构建需要已知地物的精确几何结构,故从遥感数据中自动提取地物目标的几何形态成为研究热点之一。机载激光雷达(LiDAR)因能精确获取地面高程信息,在地物目标提取方面受到了广泛关注^[1-3]。当今城市居民区内建筑物密度不断增大,并混有大量绿化树木,提取过程中建筑物之间、建筑物与树木之间的相互干扰,会使该类建筑物的提取非常困难,因此研究房屋密集居民区的建筑物提取具有重要意义^[4]。

利用LiDAR提取建筑物的研究有很多,根据提取对象主要可分为三类:建筑物三维几何结构提取、建筑物数据点提取与建筑物边缘线提取。在现今大数据量、大尺度的点云数据中,建筑物边缘以较少的数据量就能充分反映区域内建筑物的平面结构信息,因而得到了广泛应用。

目前所利用的LiDAR数据类型主要包括原始点云数据、点云数据栅格化图像^[5-7]及融合点云数据和光学影像等^[8-10]。内插后的栅格化数据存在不可避免的信息损失,无法保证原点位信息的精确性^[11-12];由于难以获取同一时相的点云数据和光学影像,且两类不同数据配准困难,后两类数据的应用范围有限。故针对原始点云数据的研究对工程实践有更重要的意义。目前,研究者已经提出了很多针对原始点云数据的建筑物边缘提取方法^[13-15],主要包括基于不规则三角网的方法^[16]、传统的基于聚类^[17]与基于滤波^[18]的方法,也有一些新型算法诸如图论^[19]与随机几何^[20-22]等。它们所面临的问题大都是无法直接提取任意形状的建筑物^[23],通常都只能以直线来拟合边缘,但很多特殊建筑物是圆形或拱形等非规则形状,故上述方法无法满足任意形状的提取要求。而基于图论的方法是针对建筑物设计的,故可以较好地得到复杂建筑物的边缘数据点,但该方法未对树木的干扰进行合适的处理,这不仅会漏提建筑物与树木交界,也会错提出树木。而随机几何方法只能提取规则形状,且实现起来过于复杂,耗时较高,目前尚无法广泛应用到实践中。

针对上述问题,本文将旋转差值核估计(RDKE)^[24-25]的基本思想应用到点云数据中。RDKE在规则离散的图像域中并不能实现,而在点云数据中不但可以实现,而且可以通过将均值与方差同时作为统计测度并调整所构建对称窗口的间距,克服绿化树木遮挡建筑物的问题,提取出具有任意形状且完整的建筑物边缘。

2 算法描述

设点云数据为P={(x_i,y_i,Z_i) | (x_i,y_i)∈F,N_i =1,2,…,n},其中,i为数据点索引,n为数据点数,(x_i,y_i)为数据点i投影在地面上的坐标,Z_i为数据点i的高程,F为点云数据在地面投影覆盖的平面域。城市居民区的主要目标包括地面、建筑物和树木。在点云数据中,建筑物点和地面点的高程起伏较小,而树木点的高程起伏较大,且通常建筑物点和树木点的高程均值比其周围地面点的大。根据上述不同地物目标的数据特点,利用RDKE实现LiDAR点云建筑物的边缘点提取。

2.1 任意方向对称窗口

为计算数据点的高程均值及描述高程起伏程度的方差等统计量,需定义任意方向的对称窗口。即需要一个对点云中所有数据点具有普适性的单位距离来衡量对称窗口的尺度。因点云数据在地面投影所覆盖的区域相对较大,且包含的数据点较多,故可认为非规则分布的数据点在此区域内近似规则分布,即每个数据点在该区域内都占有一小块正方形区域^[26]。因此,可定义近似平均点间距为

d=(xmax-xmin)(ymax-ymin)n,(1)

式中:x_min、x_max、y_min和y_max分别为点云数据中所有数据点投影在地面上的最小横坐标、最大横坐标、最小纵坐标和最大纵坐标;d为每个数据点所占有的正方形区域的边长。

沿与水平线成θ[θ∈[0,π]]的检测方向构建关于数据点i对称且中心间距为2l的两圆形对称窗口w₁(c₁,r)与w₂(c₂,r),其中,l为窗口中心到数据点i的所在位置(x_i,y_i)的距离,r为窗口半径,两窗口中心c₁与c₂的横纵坐标计算公式为

xc1=xi+lcosθ,yc1=yi+lsinθxc2=xi-lcosθ,yc2=yi-lsinθ。(2)

为使窗口内数据点能充分反映区域内高程的统计特征,需保证足够的样本容量,即窗口内数据点不能太少或为空,故应取r≥d;同时,由于检测到的边缘精细程度与窗口尺寸成反比,故r的取值也不能太大,综合考虑两方面,取r=d。图1所示为对称窗口示意图。

图 1. 某一方向的对称窗口

Fig. 1. Symmetric windows in a direction

下载图片 查看所有图片

2.2 核函数

定义在数据点i两侧的对称窗口内包含的数据点上的高斯核函数为

Kxw-xi2r,yw-yi2r=0.6171exp-12xw-xi2r2exp-12yw-yi2r2,(3)

式中(x_w,y_w)为对称窗口内数据点的位置,下标w表示对称窗口。采用(3)式可对(x_w,y_w)的高程Z进行加权,使越接近数据点i的(x_w,y_w)获得的核函数值(权重)越大,反之越小。随着数据点i逐渐接近两侧具有较大高程差的边缘点,采用(3)式会使差值绝对值呈现越来越快的非线性增长趋势,最终可使利用边缘强度提取出的边缘更靠近真实边缘。

2.3 特征点提取

定义表征不同类型地物的高程平均水平及其差异性的测度,并以此作为判据提取对应的特征点。

边缘点的两侧具有较大的高程差,而高程加权均值之差可定量描述这一特性。因此,对于某数据点i,其RDKE可表示为

M1(xi,yi,θ)=∑(xw,yw)∈w1(c1,r)Z(xw,yw)Kxw-xi2r,yw-yi2r∑(xw,yw)∈w1(c1,r)Kxw-xi2r,yw-yi2r,(4)

M2(xi,yi,θ)=∑(xw,yw)∈w2(c2,r)Z(xw,yw)Kxw-xi2r,yw-yi2r∑(xw,yw)∈w2(c2,r)Kxw-xi2r,yw-yi2r,(5)

MRDKE(xi,yi,θ)=M1(xi,yi,θ)-M2(xi,yi,θ),(6)

式中M₁(x_i,y_i,θ)和M₂(x_i,y_i,θ)分别为w₁内和w₂内样本数据点高程的加权平均值,M_RDKE(x_i,y_i,θ)为它们的差值绝对值。设θ^*代表数据点i两侧具有最大高程加权均值之差的绝对值时所对应的角度,即

θ*=argmaxθ∈[0,π)[MRDKE(xi,yi,θ)]。(7)

M_RDKE(x_i,y_i,θ^*)即为数据点i两侧具有的最大高程加权均值之差的绝对值的大小,可用于衡量该数据点周围的高程突变程度。当数据点i为非边缘点时,在任何方向上的两对称窗口中所含数据点的平均高程相近,故M_RDKE(x_i,y_i,θ^*)较小;当数据点i为边缘点时,在某一方向上两对称窗口中必然分别包含高程相差较大的两类地物点,故此时边缘区域的M_RDKE(x_i,y_i,θ^*)相较其他区域的要大很多。图2所示为在边缘点处对称窗口正好取得和没有取得M_RDKE(x_i,y_i,θ^*)时的示意图。为了更好地表述M_RDKE(x_i,y_i,θ^*)及θ^*,以图3所示的测试点云数据为例来进行说明。

图 2. 正好取得与没有取得MRDKE(xi,yi,θ*)时的对称窗口示意图

Fig. 2. Schematic of symmetric windows of MRDKE(xi,yi,θ*) obtained and not obtained

下载图片 查看所有图片

图4(a)所示为不同区域的数据点对应的M_RDKE(x_i,y_i,θ^*)及对称窗口连线方向示意图,其中对称窗口颜色越浅,其中心数据点i的M_RDKE(x_i,y_i,θ^*)越小;对称窗口颜色越深,M_RDKE(x_i,y_i,θ^*)越大,两对称窗口中心的连线与水平方向夹角为θ^*。可以看出,边缘点附近的M_RDKE(x_i,y_i,θ^*)较大,且连线方向基本垂直于建筑物和树木边缘的切线方向;此外,在树木和地面区域中,数据点i周围的高程变化较小,对应的M_RDKE(x_i,y_i,θ^*)较小,且没有明显的方向性高程突变,故连线方向基本都是杂乱无章的。综上,M_RDKE(x_i,y_i,θ^*)与θ^*可以较好地定量描述点云中具有的高程差大小和窗口连线方向。

建筑物和其周围的地面都较为平坦,而树木高

图 3. 测试数据的LiDAR点云与光学影像。(a)测试数据的LiDAR点云;(b)测试数据的光学影像

Fig. 3. LiDAR point cloud data and optical image for test data. (a) LiDAR point cloud for test data; (b) optical image for test data

下载图片 查看所有图片

图 4. 测试点云数据中不同区域的数据点对应的对称窗口连线方向及MRDKE(xi,yi,θ*)和σi(θ*)大小。(a)连线方向与MRDKE(xi,yi,θ*)大小;(b)连线方向与σi(θ*)大小

Fig. 4. Directions of symmetrical window connection corresponding to the data points in different regions in the test LiDAR point cloud data and the magnitude of MRDKE(xi, yi, θ*) and σi(θ*). (a) Connection directions and the magnitude of MRDKE(xi, yi, θ*); (b) connection directions and the magnitude of σi(θ*)

下载图片 查看所有图片

低起伏较大,因此建筑物与地面交界处和树木边缘是有区别的,而高程方差之差可定量描述这一区别。对于数据点i,在其连线方向上的对称窗口内,数据点高程的对称方差之差的绝对值为

σi(θ*)=σi(w1,θ*)-σi(w2,θ*),(8)

式中σ_i(w₁,θ^*)、σ_i(w₂,θ^*)分别为窗口w₁和w₂内数据点高程的方差。σ_i(θ^*)可用于衡量在连线方向数据点i两侧高低起伏的突变程度。

图4(b)所示为不同区域的数据点对应的对称窗口连线方向及σ_i(θ^*)大小的示意图,其中对称窗口颜色越浅表示σ_i(θ^*)越小。可以看出,建筑物与地面交界处的两侧都相对平坦,因此对应的σ_i(θ^*)较小;在树木边缘处的两侧中一侧平坦而另一侧起伏很大,故σ_i(θ^*)较大;平坦地面区域中的σ_i(θ^*)较小;而在起伏很大的树木区域中,数据点微小的高程变动就会引起相对较大的方差变化,而每个数据点的起伏又是随机的,故该区域的σ_i(θ^*)会比地面和建筑物中的要大,而比树木边缘处的要小。综上,σ_i(θ^*)可以较好地定量描述建筑物与地面交界和树木边缘之间的区别。

为提取树木点,需要准确定量地反映树木的起伏突变程度。但σ_i(θ^*)代表的是对称窗口连线方向上起伏的突变程度,而非最大的起伏突变程度。如数据点i位于起伏较大的树木区域,但连线方向上的起伏突变很小时,σ_i(θ^*)值很小。在这种情况下,若采用σ_i(θ^*)作为判据,会造成树木点的漏检,图4(b)中树木区域内的灰色对称窗口说明了这个现象。因此,对于数据点i,其任意方向对称方差之差的绝对值为

σi(θ)=σi(w1,θ)-σi(w2,θ)。(9)

设θ^#代表在数据点i处所有对称窗口连线方向中最大高程起伏突变所对应的角度,则有

θ#=argmaxθ∈[0,π)[σi(θ)]。(10)

σ_i(θ^#)可衡量数据点i周围高低起伏的最大突变程度。因方差对高程变动非常敏感,故树木区域的σ_i(θ^#)会始终比平坦区域的σ_i(θ^#)较大,最终使得树木处的每个树木点都有相对于平坦区域的较大σ_i(θ^#)值,故可提取完整的树木点。综上,三个用于表征不同类型地物的高程平均水平及其差异性的测度分别为M_RDKE(x_i, y_i, θ^*)、σ_i(θ^*)与σ_i(θ^#)。

为提取建筑物与地面的交界点,令l=2r,即当数据点i是以建筑物与地面交界为中心、宽度为2(l-r)的范围内的点时,两对称窗口内分别且始终包含的是较平坦的建筑物点与地面点,故此时M_RDKE(x_i, y_i, θ^*)较大且几乎不变,σ_i(θ^*)较小且几乎不变;当对称窗口开始同时包含建筑物点与地面点时,σ_i(θ^*)逐渐变大,而当它们的数量正好相同时,σ_i(θ^*)达到最大;当对称窗口内建筑物点和地面点的数量又逐渐变得不同时,σ_i(θ^*)逐渐变小。图5所示为当l=2r时数据点i及对称窗口穿过建筑物与地面交界的示意图,且图6(a)和图6(b)所示分别为对应的M_RDKE(x_i, y_i, θ^*)和σ_i(θ^*)的变化趋势示意图。

图 5. 数据点i和窗口的移动示意图。(a) θ*=0且l=2r时的对称窗口;(b)数据点i经过建筑物与地面交界

Fig. 5. Diagram of movement of data point i and windows. (a) Symmetric windows with θ*=0 and l=2r; (b) data point i through the boundary between building and ground

下载图片 查看所有图片

图 6. MRDKE(xi, yi, θ*)与σi(θ*)的变化趋势。(a) MRDKE(xi, yi, θ*);(b) σi(θ*)

Fig. 6. Change trend of MRDKE(xi, yi, θ*) and σi(θ*). (a) MRDKE(xi, yi, θ*); (b) σi(θ*)

下载图片 查看所有图片

图 7. 与图6中MRDKE(xi, yi, θ*)与σi(θ*)的变化趋势对应的整体点云示意图。(a)与图6(a)中宽度为2(l-r)且MRDKE(xi, yi, θ*)较大的一段红线相对应;(b1)与图6(b)中宽度为2(l-r)且σi(θ*)较小的一段红线相对应; (b2)与图6(b)中两个凸起的峰的顶端相对应

Fig. 7. Schematic of overall point cloud corresponding to the change trend of MRDKE(xi, yi, θ*) and σi(θ*) in Fig. 6. (a) Corresponding to a red line with a width of 2(l-r) and a larger MRDKE(xi, yi, θ*) in Fig. 6(a); (b1) corresponding to a red line with a width of 2(l-r) and a smaller σi(θ*) in Fig. 6(b); (b2) Corresponding to the top of the two raised peak in Fig. 6(b)

下载图片 查看所有图片

与图6中M_RDKE(x_i,y_i,θ^*)与σ_i(θ^*)的变化趋势对应的整体点云示意图如图7所示。图7(a)和图7(b1)中建筑物与地面交界两侧的白点和黑点分别与图6(a)中M_RDKE(x_i,y_i,θ^*)曲线中的X₁和X₂和图6(b)中σ_i(θ^*)曲线中的X₃和X₄对应,两位置的间距都为2(l-r)。图7(b2)表示窗口同时包含了几乎等量的建筑物点与地面点,σ_i(θ^*)在此处取极大值,建筑物与地面交界两侧的白点对应于图6(b)中两个凸起的峰的顶端。

综上,建筑物与地面交界附近的点同时具有较大的M_RDKE(x_i,y_i,θ^*)与较小的σ_i(θ^*),明显区别于其他点。故当M_RDKE(x_i,y_i,θ^*)大于阈值T₁且σ_i(θ^*)小于阈值T₂时,数据点i为建筑物与地面交界点;否则为其他点。最终可检测出以建筑物与地面交界为中心、宽度比2(l-r)稍宽的范围内的点。

建筑物边缘由建筑物与地面交界和建筑物与树木交界组成,但上述操作无法识别建筑物与树木交界。故可先将树木点提取出来,并在原点云数据上进行滤除,再提取建筑物边缘。接下来进行树木点的提取。

当l=2r且数据点i是以树木边缘为中心、2(l-r)宽度范围内的点时,两对称窗口内分别且始终包含的是树木点与地面点或建筑物点,此时σ_i(θ^#)很大且几乎不变;而当数据点i位于树木内部区域时,两窗口内分别且始终包含的是起伏同样较大的树木点,这时的σ_i(θ^#)远小于树木边缘处的σ_i(θ^#)。此时若进行阈值处理,就易将一些树木内部区域的点漏提。

因此,为提取完整的树木点,宽度调整为l=r。当数据点i靠近树木边缘时,其中一个窗口会开始同时包含树木点和地面点,此时σ_i(θ^#)开始比平坦区域的大,但其并不像l=2r且在2(l-r)宽度范围内时会全部包含树木点,故当l=r时,树木边缘处的σ_i(θ^#)没有当l=2r时树木边缘处的大。当l=r时,树木边缘处的σ_i(θ^#)与树木内部区域的σ_i(θ^#)相差不大,这会避免漏提一些树木内部区域的数据点。

图8所示为当l=r时数据点i穿过树木以及对应的σ_i(θ^#)变化趋势示意图。综上可知,阈值处理σ_i(θ^#)后,提取结果中可有效保留树木点。因此,当σ_i(θ^#)大于阈值T₃时,可判定数据点i为树木点;否则数据点i为其他点。

图 8. 数据点i穿过树木的示意图及对应的σi(θ#)的变化趋势。(a)数据点i在树木处;(b)数据点i穿过树木;(c) σi(θ#)的变化趋势

Fig. 8. Diagrammatic sketch of data point i passing through the tree and the change trend of σi(θ#).(a) Data point i at tree; (b) data point i passing through tree; (c) change trend of σi(θ#)

下载图片 查看所有图片

图 9. 数据点i穿过建筑物与地面交界的示意图及穿过该交界时σi(θ#)的变化趋势。(a) θ*=0且l=r时的对称窗口;(b)数据点i经过建筑物与地面交界的过程示意图;(c)σi(θ#)的变化趋势

Fig. 9. Diagrammatic sketch of data point i through the boundary between building and ground and the corresponding change trend to σi(θ#). (a) Symmetric windows when θ*=0 and l=r; (b) schematic of the process of data point i through the boundary between building and ground; (c) change trend of σi(θ#)

下载图片 查看所有图片

然而,当l=r且数据点i通过建筑物与地面交界时,窗口中会开始同时包含建筑物点和地面点,σ_i(θ^#)逐渐变大。故经阈值处理σ_i(θ^#)后,除了会提取出树木点,也会提取出以建筑物与地面交界为中心、宽度小于4l的范围内的建筑物与地面的交界点。l=r时,数据点i穿过建筑物与地面交界以及对应的σ_i(θ^#)变化趋势示意图如图9所示,可以看出,σ_i(θ^#)的变化趋势与l=2r时的相似,但不存在σ_i(θ^*)中2(l-r)宽度范围内的一段不变的较小值。

图 10. l=r时,不同区域的数据点对应的对称窗口连线方向及σi(θ#)大小

Fig. 10. Directions of symmetric window connection and magnitude of σi(θ#) corresponding to the data points in different regions when l=r

下载图片 查看所有图片

l=r时,不同区域的数据点对应的对称窗口连线方向及σ_i(θ^#)大小示意图如图10所示。对称窗口的颜色越浅表示它们的中心数据点i的σ_i(θ^#)越小。而在建筑物与地面交界附近的深色对称窗口表示当数据点i靠近该交界时,某一个窗口同时包含了几乎等量的建筑物点与地面点,此处的σ_i(θ^#)有极大值,即建筑物与地面交界两侧的白点所在的位置,也对应于图9(c)中两个凸起的峰的顶端。还可直观看出,建筑物与地面交界附近的σ_i(θ^#)和树

木区域的σ_i(θ^#)都较大,故阈值处理σ_i(θ^#)后提取时会将树木点和建筑物与地面交界附近的点同时提取出来。

为提取树木点,需将同时提取出的建筑物与地面的交界点删除。设利用M_RDKE(x_i, y_i, θ^*)和σ_i(θ^*)提取出的建筑物与地面交界点集为A,利用σ_i(θ^#)提取出的树木及建筑物和地面交界点集为B,用O表示将A在B中对应的点都删除后的结果。A是以建筑物与地面交界为中心、宽度比2(l-r)=2r(此时l=2r)稍宽的范围内的点集,B中的建筑物与地面交界点集是以建筑物与地面交界为中心、宽度小于4l=4r(此时l=r)的范围内的点。故在O中的建筑物与地面交界两侧会残留一些总宽度小于4r-2r=2r的少量点集,而单侧的宽度小于r。图11所示为经阈值处理后的点集范围示意图,褐色与绿色实线分别为经过T₁和T₂阈值处理后的点集范围;紫色实线为二者同时处理后的点集范围;黄色实线为经过T₃阈值处理后的点集范围;红色实线为建筑物与地面交界两侧的少量点集范围。

图 11. 阈值处理后的点集范围。(a) T1和T2;(b) T3

Fig. 11. Range of point sets after threshold procession. (a) T1 and T2; (b) T3

下载图片 查看所有图片

由于O是二值化点集,即树木点及建筑物和地面交界两侧的少量点集的高程为1,其他点的高程为0,故可通过滤波操作将建筑物与地面交界两侧的少量点集滤除:

fGxi,yi=Gxi,yi,t1≥t0G0xi,yi,t1<t0,11

Of={Of(xi,yi)|(xi,yi)∈F},(12)

式中:G(x_i,y_i)为在O中以高程为1的数据点i为圆心、2r为半径的圆形区域内的所有点的集合;t₁和t₀分别为G(x_i,y_i)中高程为1和高程为0的点数;G₀(x_i, y_i)表示G(x_i,y_i)中所有点的高程都置0后的点集;O_f为对O进行(11)式所示的滤波操作后的结果。

2.4 建筑物边缘提取

O_f中主要为树木点,且激光一般会透过树冠打在地面上,故点云中大量地面点与树木点会混杂在一起。大量地面点的高程相差很小,故将树木所在区域一定范围内的所有数据点的高程中出现次数最多的值近似为此区域的地面高程。图12所示为以某树木区域为例进行树木点滤除的过程,由图12(c)可知,5约是出现次数最多的高程值。因此,将所有点的高程都置为该高程值,所有点的高程都变为了此区域的近似地面高程(5左右)。同理,对O_f进行规则分块,对每个单位块执行上述操作,得到滤除了O_f后的结果P_h,再对P_h的每个单位块分别进行建筑物与地面交界点的提取。不对整个P_h进行提取是为防止两个相邻的块由于高程相差过大而误将它们的边界提取出来。由于已经滤除了树木点,此时的建筑物与地面的交界点即为建筑物边缘点。

图 12. 滤除树木的过程示意图。(a)树木点的三维显示;(b)树木点的二维显示;(c)高程直方图;(d)处理后的树木点二维显示

Fig. 12. Diagram of the process of filtering trees. (a) Three dimensional version of tree points; (b) two dimensional version of tree points; (c) elevation histogram; (d) two dimensional version of treated tree points

下载图片 查看所有图片

可以通过先验知识来人工选取阈值,也可通过下述方案来确定。由于建筑物边缘点和树木点等只占所有数据点的小部分,而绝大部分的数据点都位于平坦区域,故遍历点云后得到的σ_i(θ^*)中绝大部分的值都集中在较小的值附近(对应平坦区域),只有很少一部分在较大的值中均匀分布(对应建筑物边缘点和树木点)。在σ_i(θ^*)的直方图中,较小的值和较大的值之间的临界点所对应的值就是T₂。图13所示为测试数据的σ_i(θ^*)与σ_i(θ^#)的直方图。从直方图中的第一条起算,开始时值逐渐递减,到临界点时开始增大,故认为该点为T₂。同理,T₃也需根据σ_i(θ^#)的直方图确定,而T₁可根据经验选取。

图 13. 测试数据的σi(θ*)与σi(θ#)的直方图与阈值。(a) σi(θ*)的直方图;(b)图13(a)的局部放大;(c) σi(θ#)的直方图;(d)图13(c)的局部放大

Fig. 13. Histograms and thresholds of σi(θ*) and σi(θ#). (a) Histogram of σi(θ*); (b) local magnification of Fig.13(a); (c) histogram of σi(θ#); (d) local magnification of Fig.13(c)

下载图片 查看所有图片

2.5 算法流程

为尽可能实现理论上的任意方向对称窗口,以Δθ为角度间隔使θ遍历[0, π],可产生[(π-Δθ)/Δθ]+1个不同的角度θ_j=(j-1)π/{[(π-Δθ)/Δθ]+1},j=1,2,…,[(π-Δθ)/Δθ]+1,其中j为角度索引。将[(π-Δθ)/Δθ]+1设为k,即可产生k个不同角度对应的对称窗口。所提算法的流程可总结如下。

S1 计算点云数据的近似点间距d,确定窗口尺寸,在数据点i上构建k个不同角度的对称窗口。

S2 通过(4)~(8)式得到M_RDKE(x_i,y_i,θ^*)和σ_i(θ^*),并对所有点云数据都进行该操作。

S3 将每个数据点的M_RDKE(x_i,y_i,θ^*)和σ_i(θ^*)通过T₁和T₂阈值处理后得到建筑物与地面交界A。

S4 将每个数据点的σ_i(θ^#)通过T₃阈值处理后得到树木及建筑物和地面交界B。

S5 在B中删除A中对应位置的建筑物与地面交界后得到O,并将建筑物与地面交界两侧的少量点集滤除后得到O_f。

S6 对O_f进行规则分块,滤除O_f后得到结果P_h,再对P_h进行规则分块并逐块提取出建筑物边缘。

设A中建筑物与地面交界点的总数为p,O_f中检测出的树木点总数为q。由S1和S2可知,d的计算与其他变量无关,故其复杂度为O(1)。M_RDKE(x_i,y_i,θ^*)和σ_i(θ^*)需要遍历k个不同角度的对称窗口后才能求得,同时还要遍历点云数据,即再循环n次,故S1和S2的总复杂度为O(1)+O(kn)。S3为将阈值操作遍历点云,其复杂度为O(n)。S4的第一步需要计算所有数据点的σ_i(θ^#),其复杂度为O(kn),第二步的复杂度为O(n),故S4的复杂度为O(kn)+O(n);S5需要删除p个点,还要遍历点云来滤波,故S5的复杂度为O(p)+O(n);同理可得,S6的复杂度为O(q)+O(n),其中q为S6需要删除的点个数。

综上,所提算法的复杂度为O(1)+O(kn)+O(n)+O(kn)+O(n)+O(p)+O(n)+O(q)+O(n),简化后为O(2kn+4n+p+q+1),只保留最高阶项且去除常系数,故最终复杂度为O(kn)。当k取一个固定的值时,可将其视为常数,此时所提算法的复杂度为O(n)。

3 实验结果与讨论

应用所提算法对测试数据进行提取,图14所示为建筑物与地面交界A的提取过程。在图14(a)中,所有边缘点的边缘强度都较大,且图14(b)表明所有边缘点都以最大高程差突变的方向被检测到,从而使任意方向的边缘强度都一致而便于提取;而在图14(c)中,只有建筑物与地面交界处的边缘强度较小,故将图14(a)和(c)结合后可以较好地分离出图14(e)中的建筑物与地面交界。图15所示为树木及建筑物和地面交界B的提取过程,可以看出,通过σ_i(θ^#)可以尽可能地将所有树木点的强度提高到相对于平坦区域的较大值,以便于提取。

图 14. 建筑物与地面交界的提取过程。(a) MRDKE(xi, yi, θ*)量级分布图;(b) MRDKE(xi, yi, θ*)及θ*所构成的矢量场;(c) σi(θ*)量级分布图;(d)σi(θ*)及θ*所构成的矢量场;(e)提取出的建筑物与地面交界

Fig. 14. Extraction process of the boundaries between buildings and ground. (a) Magnitude distribution of MRDKE(xi, yi, θ*); (b) vector field consisting of MRDKE(xi, yi, θ*) and θ*; (c) magnitude distribution of σi(θ*); (d) vector field consisting of σi(θ*) and θ*; (e) the extracted boundaries between buildings and ground

下载图片 查看所有图片

图 15. 树木及建筑物和地面交界的提取过程示意图。(a) σi(θ#)量级分布图;(b) σi(θ#)及θ#所构成的矢量场;(c)提取出的树木及建筑物和地面交界

Fig. 15. Schematic of the extraction process of trees and the boundaries between buildings and ground. (a) Magnitude distribution of σi(θ#); (b) vector field consisting of σi(θ#) and θ#; (c) trees and buildings and the boundaries between ground

下载图片 查看所有图片

图16所示为最终的建筑物边缘提取过程。由图16(a)和(b)可知,滤除建筑物与地面交界两侧的少量点集后,可以较好地保留树木点;由图16(c)和(d)可知,将O_f进行规则分块可以较好地将树木点的高程还原为其所在地面区域的高程,进而起到滤除树木的作用。

为判断所提算法的准确性,以标准边缘点在算法提取出的边缘点中所占的比率作为定量评价标准,将以人为确定的边缘线为中心、半径为d的范围内包含的所有数据点作为标准边缘点。表1所示

图 16. 建筑物边缘的提取过程。(a)点集O;(b)点集Of;(c) Of的规则块划分;(d)点集Ph;(e)建筑物边缘

Fig. 16. Extraction process of building edges. (a) Point set O; (b) point set Of; (c) regular block partition of Of; (d) point set Ph; (e) building edges

下载图片 查看所有图片

为定量评价结果,可以看出,所提算法的平均提取率可达到80%左右。

为进一步验证所提算法的实用性和有效性,分别对图17所示的4个大尺度数据进行实验。图18所示为所提算法应用到4个实验数据后的A、B、O_f与P_h的结果。为了证明所提算法相对于其他方法的有效性,将通过Terrasolid软件的基于滤波的方法、文献[ 19]中基于图论的方法与文献[ 20]中基于随机几何的方法作为对比算法,将所提算法和对比算法的提取结果进行对比。图19所示为所提算法的提取结果,图20所示为三个对比算法的提取结果。

图 17. 居民区点云数据及其对应的光学影像。(a)~(d)点云数据;(a1)~(d1)与点云数据对应的光学影像

Fig. 17. Point cloud data of residential areas and corresponding optical images. (a)-(d) Point cloud data; (a1)-(d1) optical images corresponding to the point cloud data

下载图片 查看所有图片

为了进行定性评价,分别将所提算法与对比算法的提取结果叠加到原始点云数据上。图21所示为所提算法的叠加结果,图22所示为对比算法的叠加结果。

由图22可知,Terrasolid软件仅以几何多边形来拟合其经过滤波后的建筑物点,这就造成了相隔非常近的建筑物区分不开,且无法很好地拟合非规则建筑物,难以得到理想的提取结果;文献[ 19]中的方法并没有针对树木准确定义相应的高程结构模型,从而很容易误追踪出与建筑物相粘连的树木边缘,这大大降低了提取的准确性;文献[ 20]中的方法仅以矩形标识来拟合建筑物,这直接导致矩形无法完全“套”住不规则建筑物。更严重的是,该方法求解矩形标识的角度不准确,矩形经常不能按照建筑物的方向来拟合,使准确性降低。且当建筑物较大时,其只能以较小尺寸的矩形来拼凑着拟合大型建

图 18. 第一个到第4个实验数据对应的A、B、Of、Ph。(a1)~(d1) A;(a2)~(d2) B;(a3)~(d3) Of;(a4)~(d4) Ph

Fig. 18. A、B、Of、Ph for first to fourth experimental data. (a1)-(d1) A; (a2)-(d2) B; (a3)-(d3) Of; (a4)-(d4) Ph

下载图片 查看所有图片

图 19. 所提算法的提取结果。(a)数据一;(b)数据二;(c)数据三;(d)数据四

Fig. 19. Extraction results based on proposed algorithm. (a) Data 1; (b) data 2; (c) data 3; (d) data 4

下载图片 查看所有图片

图 20. 对比算法的提取结果。(a1)~(d1) Terrasolid软件;(a2)~(d2)文献[ 19]中的方法;(a3)~(d3)文献[ 20]中的方法

Fig. 20. Extraction results based on compared methods. (a1)-(d1) Terrasolid software; (a2)-(d2) the method in Ref. [19]; (a3)-(d3) the method in Ref. [20]

下载图片 查看所有图片

图 21. 所提算法的叠加结果。(a)数据一;(b)数据二;(c)数据三;(d)数据四

Fig. 21. Superposition results of proposed algorithm. (a) Data 1; (b) data 2; (c) data 3; (d) data 4

下载图片 查看所有图片

图 22. 对比算法的叠加结果。(a1)~(d1)Terrasolid软件;(a2)~(d2)文献[ 19]中的方法;(a3)~(d3)文献[ 20]中的方法

Fig. 22. Superposition results of compared methods. (a1)-(d1) Terrasolid software; (a2)-(d2) the method in Ref. [19]; (a3)-(d3) the method in Ref. [20]

下载图片 查看所有图片

筑物,导致很多“漏洞”未被拟合,不适用于工程实践。

而由图21可知,所提算法不但能够区分开相隔非常近的建筑物,而且可以提取出不规则的边缘。不仅如此,所提算法可以较好地滤除掉树木点,克服了传统方法错提与漏提建筑物的问题,满足了工程应用上的实际需求。

4 结论

为克服建筑物边缘提取方法中无法检测任意形状建筑物边缘,且受树木遮挡使得检测不完整的问题,通过在点云数据中构建任意方向的对称窗口,将RDKE引入到点云数据中。通过将加权均值与方差同时作为统计测度,调整窗口间距,提取出树木点,并利用激光穿透树冠的特性,对树木点进行规则分块,完整地还原了点云中树木所处位置的地面高程信息,滤除了树木点,解决了树木阻挡建筑物导致的无法利用RDKE直接提取完整建筑物边缘的问题,提取出了任意几何形状的建筑物边缘。但所提算法提取的是建筑物边缘处的点集,并不是完全精确的单个边缘点,后续会进一步改进该问题。